Working ROCI analysis

#RWarngle data
med_hist <- med_hist %>% dplyr::select(subj_id, heart_disease___2, copd) %>% rename(Subject_ID = subj_id) %>% mutate(chf = as.factor(heart_disease___2), copd = as.factor(copd)) %>% dplyr::select(-heart_disease___2)

outcomes <- outcomes %>% mutate(non_white = as.factor(ifelse(Race == "White, non-Hispanic origin", 0, 1)), hosp_death = as.factor(ifelse(DC_coded == 0 | DC_coded == 1, 1, 0)), Immsupp = as.factor(Immsupp), Male = as.factor(Male)) %>% left_join(med_hist, by = "Subject_ID")

imgs <- imgs %>% dplyr::select(-day_3_cxr, -day_5_cxr) %>% filter(admit_cxr == 1) %>% left_join(outcomes, by = "Subject_ID") %>% filter(Final_Dx != "Control") 

#Numerics to factor
imgs$mort28d <- as.factor(imgs$mort28d)
imgs$mort60d <- as.factor(imgs$mort60d)
rm(outcomes)
rm(med_hist)

refdate <- as.Date("2018-09-15")
#Create variables with survival (in days) and a binary indicating censored status
imgs <- mutate(imgs, survival = ifelse(is.na(date_death) == TRUE, as.numeric(difftime(refdate, ICU_admit, units='days')), as.numeric(difftime(date_death, ICU_admit, units='days')))) %>% mutate(censor = ifelse(is.na(date_death) == TRUE, 0, 1))

In-hospital Mortality

Logistic Regression

#In-hospital mortality
hosp_mort <- glm(hosp_death ~ Total+Male+Age+Immsupp+APACHE+non_white, family=binomial(link='logit'), data=imgs)
summary(hosp_mort)

## 
## Call:
## glm(formula = hosp_death ~ Total + Male + Age + Immsupp + APACHE + 
##     non_white, family = binomial(link = "logit"), data = imgs)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6980  -0.7684  -0.5036   0.9073   2.3924  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -4.084395   0.594178  -6.874 6.24e-12 ***
## Total        0.092477   0.024730   3.739 0.000184 ***
## Male1        0.053696   0.209066   0.257 0.797303    
## Age          0.006444   0.006662   0.967 0.333365    
## Immsupp1     1.269744   0.213219   5.955 2.60e-09 ***
## APACHE       0.060922   0.013030   4.676 2.93e-06 ***
## non_white1  -0.551567   0.291906  -1.890 0.058820 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 668.22  on 559  degrees of freedom
## Residual deviance: 563.83  on 553  degrees of freedom
## AIC: 577.83
## 
## Number of Fisher Scoring iterations: 4

cbind(OR = exp(coef(hosp_mort)),
      exp(confint(hosp_mort)),
      p = coef(summary(hosp_mort))[,4])

##                     OR       2.5 %     97.5 %            p
## (Intercept) 0.01683332 0.005062654 0.05219008 6.241538e-12
## Total       1.09688830 1.045465785 1.15209052 1.844135e-04
## Male1       1.05516410 0.700653326 1.59203449 7.973030e-01
## Age         1.00646525 0.993478299 1.01981579 3.333653e-01
## Immsupp1    3.55994085 2.351709648 5.43100824 2.599010e-09
## APACHE      1.06281649 1.036392348 1.09082005 2.929309e-06
## non_white1  0.57604650 0.318031588 1.00427699 5.882011e-02

#In-hospital mortality, entire cohort
hosp_mort_phx <- glm(hosp_death ~ Total+Male+Age+Immsupp+APACHE+non_white+copd+chf, family=binomial(link='logit'), data=imgs)
summary(hosp_mort)

## 
## Call:
## glm(formula = hosp_death ~ Total + Male + Age + Immsupp + APACHE + 
##     non_white, family = binomial(link = "logit"), data = imgs)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6980  -0.7684  -0.5036   0.9073   2.3924  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -4.084395   0.594178  -6.874 6.24e-12 ***
## Total        0.092477   0.024730   3.739 0.000184 ***
## Male1        0.053696   0.209066   0.257 0.797303    
## Age          0.006444   0.006662   0.967 0.333365    
## Immsupp1     1.269744   0.213219   5.955 2.60e-09 ***
## APACHE       0.060922   0.013030   4.676 2.93e-06 ***
## non_white1  -0.551567   0.291906  -1.890 0.058820 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 668.22  on 559  degrees of freedom
## Residual deviance: 563.83  on 553  degrees of freedom
## AIC: 577.83
## 
## Number of Fisher Scoring iterations: 4

cbind(OR = exp(coef(hosp_mort_phx)),
      exp(confint(hosp_mort_phx)),
      p = coef(summary(hosp_mort_phx))[,4])

##                     OR       2.5 %     97.5 %            p
## (Intercept) 0.01588295 0.004742023 0.04951704 4.085487e-12
## Total       1.10425448 1.051822966 1.16068666 7.672017e-05
## Male1       1.03405987 0.684772429 1.56386419 8.734939e-01
## Age         1.00820910 0.994906152 1.02186938 2.298141e-01
## Immsupp1    3.54015781 2.331276815 5.41868157 4.006920e-09
## APACHE      1.06259782 1.036062340 1.09073793 3.517687e-06
## non_white1  0.59180355 0.325476766 1.03631651 7.437630e-02
## copd1       0.88654597 0.485450758 1.58012353 6.880652e-01
## chf1        0.47009908 0.196615699 1.02947211 7.118681e-02

imgs_non <- imgs %>% filter(Final_Dx != "ARDS" & Final_Dx != "Sepsis/ARDS") %>% mutate(quartile = as.factor(ntile(Total,4)))

hosp_mort_non <- glm(hosp_death ~ Total+Male+Age+Immsupp+APACHE+non_white+copd+chf, family=binomial(link='logit'), data=imgs_non)
summary(hosp_mort_non)

## 
## Call:
## glm(formula = hosp_death ~ Total + Male + Age + Immsupp + APACHE + 
##     non_white + copd + chf, family = binomial(link = "logit"), 
##     data = imgs_non)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.4545  -0.7166  -0.4920  -0.2853   2.4353  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.539439   0.679278  -5.211 1.88e-07 ***
## Total        0.079191   0.030246   2.618 0.008840 ** 
## Male1       -0.101988   0.251225  -0.406 0.684770    
## Age          0.006954   0.008245   0.843 0.399012    
## Immsupp1     0.897279   0.258076   3.477 0.000507 ***
## APACHE       0.049636   0.016075   3.088 0.002017 ** 
## non_white1  -0.884079   0.373005  -2.370 0.017781 *  
## copd1       -0.093300   0.343324  -0.272 0.785810    
## chf1        -0.920783   0.524107  -1.757 0.078942 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 451.96  on 435  degrees of freedom
## Residual deviance: 401.46  on 427  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: 419.46
## 
## Number of Fisher Scoring iterations: 5

cbind(OR = exp(coef(hosp_mort_non)),
      exp(confint(hosp_mort_non)),
      p = coef(summary(hosp_mort_non))[,4])

##                     OR       2.5 %    97.5 %            p
## (Intercept) 0.02902962 0.007300731 0.1053285 1.882420e-07
## Total       1.08241058 1.020403350 1.1491916 8.839609e-03
## Male1       0.90304061 0.551500761 1.4798499 6.847704e-01
## Age         1.00697821 0.990945005 1.0235943 3.990123e-01
## Immsupp1    2.45292020 1.480844976 4.0814016 5.074339e-04
## APACHE      1.05088820 1.018520687 1.0849837 2.016853e-03
## non_white1  0.41309452 0.188756381 0.8261595 1.778076e-02
## copd1       0.91091977 0.453176072 1.7541579 7.858096e-01
## chf1        0.39820728 0.127265437 1.0316932 7.894192e-02

imgs_ards <- imgs %>% filter(Final_Dx == "ARDS" | Final_Dx == "Sepsis/ARDS") %>% mutate(quartile = as.factor(ntile(Total,4)))

hosp_mort_ards <- glm(hosp_death ~ Total+Male+Age+Immsupp+APACHE+non_white+copd+chf, family=binomial(link='logit'), data=imgs_ards)
summary(hosp_mort_ards)

## 
## Call:
## glm(formula = hosp_death ~ Total + Male + Age + Immsupp + APACHE + 
##     non_white + copd + chf, family = binomial(link = "logit"), 
##     data = imgs_ards)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.2178  -0.7867   0.5077   0.8281   1.7691  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -4.84621    1.44806  -3.347 0.000818 ***
## Total        0.07842    0.06488   1.209 0.226815    
## Male1        0.46943    0.45559   1.030 0.302843    
## Age          0.02932    0.01466   2.000 0.045461 *  
## Immsupp1     2.14922    0.45618   4.711 2.46e-06 ***
## APACHE       0.04088    0.02895   1.412 0.157902    
## non_white1   0.26525    0.67167   0.395 0.692910    
## copd1        0.29880    0.69194   0.432 0.665867    
## chf1         0.35172    0.97860   0.359 0.719291    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 169.86  on 122  degrees of freedom
## Residual deviance: 133.67  on 114  degrees of freedom
## AIC: 151.67
## 
## Number of Fisher Scoring iterations: 4

cbind(OR = exp(coef(hosp_mort_ards)),
      exp(confint(hosp_mort_ards)),
      p = coef(summary(hosp_mort_ards))[,4])

##                     OR        2.5 %     97.5 %            p
## (Intercept) 0.00785811 0.0003813585  0.1156035 8.178336e-04
## Total       1.08157437 0.9541773151  1.2329778 2.268146e-01
## Male1       1.59907464 0.6563719415  3.9592858 3.028427e-01
## Age         1.02975509 1.0014229410  1.0611299 4.546108e-02
## Immsupp1    8.57816955 3.6322592529 21.9605301 2.460273e-06
## APACHE      1.04172710 0.9854780339  1.1050301 1.579016e-01
## non_white1  1.30375473 0.3516770076  5.0506735 6.929104e-01
## copd1       1.34823571 0.3458602420  5.3895929 6.658666e-01
## chf1        1.42150455 0.2154587065 11.0495013 7.192914e-01

28-day and 60-day mortality

Logistic Regression

one_month <- glm(mort28d ~ Total+Male+Age+Immsupp+APACHE+non_white, family=binomial(link='logit'), data=imgs)
summary(one_month)

## 
## Call:
## glm(formula = mort28d ~ Total + Male + Age + Immsupp + APACHE + 
##     non_white, family = binomial(link = "logit"), data = imgs)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7218  -0.7811  -0.4991   0.9232   2.4334  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -4.266691   0.601511  -7.093 1.31e-12 ***
## Total        0.091542   0.024771   3.696 0.000219 ***
## Male1        0.217256   0.209525   1.037 0.299783    
## Age          0.011328   0.006706   1.689 0.091206 .  
## Immsupp1     1.161740   0.213160   5.450 5.03e-08 ***
## APACHE       0.056919   0.012974   4.387 1.15e-05 ***
## non_white1  -0.825997   0.306350  -2.696 0.007012 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 666.36  on 559  degrees of freedom
## Residual deviance: 564.09  on 553  degrees of freedom
## AIC: 578.09
## 
## Number of Fisher Scoring iterations: 5

cbind(OR = exp(coef(one_month)),
      exp(confint(one_month)),
      p = coef(summary(one_month))[,4])

##                     OR       2.5 %     97.5 %            p
## (Intercept) 0.01402812 0.004152741 0.04405301 1.309565e-12
## Total       1.09586263 1.044406767 1.15110288 2.194154e-04
## Male1       1.24266251 0.825248052 1.87846009 2.997831e-01
## Age         1.01139207 0.998293195 1.02493661 9.120584e-02
## Immsupp1    3.19548766 2.110319502 4.87225504 5.034364e-08
## APACHE      1.05856976 1.032332806 1.08630840 1.148225e-05
## non_white1  0.43779840 0.233268832 0.78054470 7.012449e-03

two_month <- glm(mort60d ~ Total+Male+Age+Immsupp+APACHE+non_white, family=binomial(link='logit'), data=imgs)
summary(two_month)

## 
## Call:
## glm(formula = mort60d ~ Total + Male + Age + Immsupp + APACHE + 
##     non_white, family = binomial(link = "logit"), data = imgs)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.8530  -0.8407  -0.5037   0.9667   2.4456  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.78109    0.56629  -6.677 2.44e-11 ***
## Total        0.10888    0.02392   4.551 5.33e-06 ***
## Male1        0.12603    0.20222   0.623 0.533145    
## Age          0.01237    0.00644   1.921 0.054753 .  
## Immsupp1     1.21045    0.20663   5.858 4.68e-09 ***
## APACHE       0.04404    0.01247   3.533 0.000411 ***
## non_white1  -0.97089    0.29281  -3.316 0.000914 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 710.56  on 559  degrees of freedom
## Residual deviance: 595.10  on 553  degrees of freedom
## AIC: 609.1
## 
## Number of Fisher Scoring iterations: 4

cbind(OR = exp(coef(two_month)),
      exp(confint(two_month)),
      p = coef(summary(two_month))[,4])

##                     OR       2.5 %     97.5 %            p
## (Intercept) 0.02279789 0.007272274 0.06717893 2.438953e-11
## Total       1.11502549 1.064481362 1.16931188 5.333626e-06
## Male1       1.13431354 0.763559242 1.68874276 5.331450e-01
## Age         1.01244742 0.999846131 1.02545694 5.475289e-02
## Immsupp1    3.35499085 2.244146627 5.04950193 4.681115e-09
## APACHE      1.04502604 1.020052562 1.07124060 4.106833e-04
## non_white1  0.37874625 0.208044832 0.65943028 9.140911e-04

imgs %>% filter(mort60d == 1) %>% summarize(count = length(Total))

## # A tibble: 1 x 1
##   count
##   <int>
## 1   185

Overall Mortality

Cox regression of entire cohort (not in quartiles)

#Cox regression model
cox_all <- coxph(Surv(survival, censor) ~ Male+Age+non_white+APACHE+Immsupp+Total, data=imgs)
summary(cox_all)

## Call:
## coxph(formula = Surv(survival, censor) ~ Male + Age + non_white + 
##     APACHE + Immsupp + Total, data = imgs)
## 
##   n= 560, number of events= 369 
## 
##                 coef exp(coef)  se(coef)      z Pr(>|z|)    
## Male1       0.055338  1.056898  0.105766  0.523  0.60082    
## Age         0.017286  1.017436  0.003438  5.027 4.97e-07 ***
## non_white1 -0.344413  0.708636  0.141348 -2.437  0.01482 *  
## APACHE      0.026503  1.026857  0.006851  3.868  0.00011 ***
## Immsupp1    0.965437  2.625934  0.110705  8.721  < 2e-16 ***
## Total       0.053079  1.054513  0.012027  4.413 1.02e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##            exp(coef) exp(-coef) lower .95 upper .95
## Male1         1.0569     0.9462    0.8590    1.3004
## Age           1.0174     0.9829    1.0106    1.0243
## non_white1    0.7086     1.4112    0.5372    0.9348
## APACHE        1.0269     0.9738    1.0132    1.0407
## Immsupp1      2.6259     0.3808    2.1137    3.2622
## Total         1.0545     0.9483    1.0299    1.0797
## 
## Concordance= 0.695  (se = 0.016 )
## Rsquare= 0.255   (max possible= 1 )
## Likelihood ratio test= 165.1  on 6 df,   p=<2e-16
## Wald test            = 158.3  on 6 df,   p=<2e-16
## Score (logrank) test = 167.7  on 6 df,   p=<2e-16

cox_all_phx <- coxph(Surv(survival, censor) ~ Male+Age+non_white+APACHE+Immsupp+Total+copd+chf, data=imgs)
summary(cox_all_phx)

## Call:
## coxph(formula = Surv(survival, censor) ~ Male + Age + non_white + 
##     APACHE + Immsupp + Total + copd + chf, data = imgs)
## 
##   n= 559, number of events= 369 
##    (1 observation deleted due to missingness)
## 
##                 coef exp(coef)  se(coef)      z Pr(>|z|)    
## Male1       0.056936  1.058588  0.106170  0.536 0.591771    
## Age         0.017627  1.017783  0.003471  5.078 3.81e-07 ***
## non_white1 -0.331374  0.717937  0.141693 -2.339 0.019353 *  
## APACHE      0.026032  1.026373  0.006864  3.792 0.000149 ***
## Immsupp1    0.960629  2.613340  0.111177  8.641  < 2e-16 ***
## Total       0.057090  1.058751  0.012197  4.681 2.86e-06 ***
## copd1       0.038081  1.038815  0.145341  0.262 0.793311    
## chf1       -0.311076  0.732658  0.201684 -1.542 0.122978    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##            exp(coef) exp(-coef) lower .95 upper .95
## Male1         1.0586     0.9447    0.8597    1.3035
## Age           1.0178     0.9825    1.0109    1.0247
## non_white1    0.7179     1.3929    0.5438    0.9478
## APACHE        1.0264     0.9743    1.0127    1.0403
## Immsupp1      2.6133     0.3827    2.1017    3.2496
## Total         1.0588     0.9445    1.0337    1.0844
## copd1         1.0388     0.9626    0.7813    1.3812
## chf1          0.7327     1.3649    0.4934    1.0879
## 
## Concordance= 0.696  (se = 0.016 )
## Rsquare= 0.259   (max possible= 1 )
## Likelihood ratio test= 167.9  on 8 df,   p=<2e-16
## Wald test            = 160.4  on 8 df,   p=<2e-16
## Score (logrank) test = 169.9  on 8 df,   p=<2e-16

cox_ards <- coxph(Surv(survival, censor) ~ Male+Age+non_white+APACHE+Immsupp+Total+copd+chf, data=imgs_ards)
summary(cox_ards)

## Call:
## coxph(formula = Surv(survival, censor) ~ Male + Age + non_white + 
##     APACHE + Immsupp + Total + copd + chf, data = imgs_ards)
## 
##   n= 123, number of events= 89 
## 
##                 coef exp(coef)  se(coef)      z Pr(>|z|)    
## Male1       0.022694  1.022954  0.231944  0.098   0.9221    
## Age         0.016957  1.017102  0.007723  2.196   0.0281 *  
## non_white1 -0.408296  0.664782  0.350104 -1.166   0.2435    
## APACHE      0.034026  1.034611  0.015078  2.257   0.0240 *  
## Immsupp1    1.504820  4.503344  0.249986  6.020 1.75e-09 ***
## Total      -0.006605  0.993416  0.032660 -0.202   0.8397    
## copd1       0.299967  1.349814  0.354405  0.846   0.3973    
## chf1        0.114992  1.121865  0.490740  0.234   0.8147    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##            exp(coef) exp(-coef) lower .95 upper .95
## Male1         1.0230     0.9776    0.6493     1.612
## Age           1.0171     0.9832    1.0018     1.033
## non_white1    0.6648     1.5043    0.3347     1.320
## APACHE        1.0346     0.9665    1.0045     1.066
## Immsupp1      4.5033     0.2221    2.7590     7.351
## Total         0.9934     1.0066    0.9318     1.059
## copd1         1.3498     0.7408    0.6739     2.704
## chf1          1.1219     0.8914    0.4288     2.935
## 
## Concordance= 0.713  (se = 0.033 )
## Rsquare= 0.357   (max possible= 0.998 )
## Likelihood ratio test= 54.41  on 8 df,   p=6e-09
## Wald test            = 48.82  on 8 df,   p=7e-08
## Score (logrank) test = 54.79  on 8 df,   p=5e-09

cox_ards_qt <- coxph(Surv(survival, censor) ~ Male+Age+non_white+APACHE+Immsupp+quartile+copd+chf, data=imgs_ards)
summary(cox_ards_qt)

## Call:
## coxph(formula = Surv(survival, censor) ~ Male + Age + non_white + 
##     APACHE + Immsupp + quartile + copd + chf, data = imgs_ards)
## 
##   n= 123, number of events= 89 
## 
##                 coef exp(coef)  se(coef)      z Pr(>|z|)    
## Male1       0.019116  1.019300  0.233130  0.082   0.9346    
## Age         0.017232  1.017382  0.007756  2.222   0.0263 *  
## non_white1 -0.407777  0.665127  0.349985 -1.165   0.2440    
## APACHE      0.035300  1.035930  0.016508  2.138   0.0325 *  
## Immsupp1    1.510938  4.530978  0.253880  5.951 2.66e-09 ***
## quartile2  -0.103852  0.901358  0.325526 -0.319   0.7497    
## quartile3  -0.076088  0.926735  0.323918 -0.235   0.8143    
## quartile4  -0.119369  0.887481  0.321491 -0.371   0.7104    
## copd1       0.294323  1.342217  0.372813  0.789   0.4298    
## chf1        0.101239  1.106541  0.493868  0.205   0.8376    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##            exp(coef) exp(-coef) lower .95 upper .95
## Male1         1.0193     0.9811    0.6454     1.610
## Age           1.0174     0.9829    1.0020     1.033
## non_white1    0.6651     1.5035    0.3350     1.321
## APACHE        1.0359     0.9653    1.0029     1.070
## Immsupp1      4.5310     0.2207    2.7548     7.452
## quartile2     0.9014     1.1094    0.4762     1.706
## quartile3     0.9267     1.0791    0.4912     1.749
## quartile4     0.8875     1.1268    0.4726     1.667
## copd1         1.3422     0.7450    0.6464     2.787
## chf1          1.1065     0.9037    0.4203     2.913
## 
## Concordance= 0.713  (se = 0.033 )
## Rsquare= 0.358   (max possible= 0.998 )
## Likelihood ratio test= 54.53  on 10 df,   p=4e-08
## Wald test            = 48.67  on 10 df,   p=5e-07
## Score (logrank) test = 54.85  on 10 df,   p=3e-08

schoenfeld_res <- cox.zph(cox_all)
plot(schoenfeld_res)

Cox regression of entire cohort by quartile

#create quartiles
imgs <- imgs %>% mutate(quartile = ntile(Total,4)) 
imgs$quartile <- as.factor(imgs$quartile)

#Cox regressions
cox_qt <- survfit(Surv(survival, censor) ~ quartile, data=imgs) #for graph
cox_qt

## Call: survfit(formula = Surv(survival, censor) ~ quartile, data = imgs)
## 
##              n events median 0.95LCL 0.95UCL
## quartile=1 140     64   3117    1802      NA
## quartile=2 140    103    206     136     526
## quartile=3 140     97    296     104     712
## quartile=4 140    105    102      38     244

#cox_quart <- coxph(Surv(survival, censor) ~ Age+Male+Immsupp+non_white+APACHE+quartile, data=imgs) #for betas
#summary(cox_quart)

cox_ards_q <- coxph(Surv(survival, censor) ~ Male+Age+non_white+APACHE+Immsupp+quartile, data=imgs_ards)
summary(cox_ards_q)

## Call:
## coxph(formula = Surv(survival, censor) ~ Male + Age + non_white + 
##     APACHE + Immsupp + quartile, data = imgs_ards)
## 
##   n= 123, number of events= 89 
## 
##                 coef exp(coef)  se(coef)      z Pr(>|z|)    
## Male1      -0.006831  0.993193  0.221369 -0.031   0.9754    
## Age         0.017974  1.018136  0.007631  2.356   0.0185 *  
## non_white1 -0.410447  0.663354  0.347343 -1.182   0.2373    
## APACHE      0.035740  1.036387  0.016526  2.163   0.0306 *  
## Immsupp1    1.478997  4.388542  0.249134  5.937 2.91e-09 ***
## quartile2  -0.119471  0.887390  0.323193 -0.370   0.7116    
## quartile3  -0.010618  0.989438  0.309665 -0.034   0.9726    
## quartile4  -0.111176  0.894781  0.318036 -0.350   0.7267    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##            exp(coef) exp(-coef) lower .95 upper .95
## Male1         0.9932     1.0069    0.6436     1.533
## Age           1.0181     0.9822    1.0030     1.033
## non_white1    0.6634     1.5075    0.3358     1.310
## APACHE        1.0364     0.9649    1.0034     1.071
## Immsupp1      4.3885     0.2279    2.6931     7.151
## quartile2     0.8874     1.1269    0.4710     1.672
## quartile3     0.9894     1.0107    0.5393     1.815
## quartile4     0.8948     1.1176    0.4797     1.669
## 
## Concordance= 0.712  (se = 0.033 )
## Rsquare= 0.355   (max possible= 0.998 )
## Likelihood ratio test= 53.91  on 8 df,   p=7e-09
## Wald test            = 48.38  on 8 df,   p=8e-08
## Score (logrank) test = 54.41  on 8 df,   p=6e-09

cox_ards_qp <- survfit(Surv(survival, censor) ~ quartile, data=imgs_ards)

plot_ards <- ggsurvplot(cox_ards_qp, data=imgs_ards, pval=TRUE, pval.size=3,
           ggtheme = theme_minimal(), legend.title="Quartile", legend.labs=c("Q1","Q2","Q3","Q4"),
           risk.table=TRUE, tables.height = 0.2, tables.theme = theme_cleantable(),
           xlim=c(0,3650), break.time.by=365.25, xscale="d_m", ylab="",
           xlab="Months", font.legend=8, title="Entire cohort", font.title=10,
           legend=c(1.1,0.6), risk.table.y.text=FALSE,
           risk.table.fontsize=3, risk.table.title="") 
plot_ards$plot <- plot_ards$plot + theme(plot.title = element_text(hjust = 0.5))
plot_ards

#Adding copd/chf did not explode the std errors. 

contrasts(imgs$quartile) <- contr.treatment(4, base=4)
cox_qt_base4 <- coxph(Surv(survival, censor) ~ Age+Male+Immsupp+non_white+APACHE+copd+chf+quartile, data=imgs)
summary(cox_qt_base4)

## Call:
## coxph(formula = Surv(survival, censor) ~ Age + Male + Immsupp + 
##     non_white + APACHE + copd + chf + quartile, data = imgs)
## 
##   n= 559, number of events= 369 
##    (1 observation deleted due to missingness)
## 
##                 coef exp(coef)  se(coef)      z Pr(>|z|)    
## Age         0.015858  1.015984  0.003519  4.506 6.60e-06 ***
## Male1       0.054028  1.055514  0.106463  0.507 0.611817    
## Immsupp1    0.995257  2.705419  0.111847  8.898  < 2e-16 ***
## non_white1 -0.347144  0.706703  0.141602 -2.452 0.014224 *  
## APACHE      0.025299  1.025622  0.006914  3.659 0.000253 ***
## copd1       0.052707  1.054120  0.147045  0.358 0.720015    
## chf1       -0.261249  0.770089  0.201493 -1.297 0.194780    
## quartile1  -0.781829  0.457568  0.163799 -4.773 1.81e-06 ***
## quartile2  -0.059434  0.942298  0.142399 -0.417 0.676403    
## quartile3  -0.105750  0.899650  0.143789 -0.735 0.462063    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##            exp(coef) exp(-coef) lower .95 upper .95
## Age           1.0160     0.9843    1.0090    1.0230
## Male1         1.0555     0.9474    0.8567    1.3004
## Immsupp1      2.7054     0.3696    2.1729    3.3685
## non_white1    0.7067     1.4150    0.5354    0.9328
## APACHE        1.0256     0.9750    1.0118    1.0396
## copd1         1.0541     0.9487    0.7902    1.4062
## chf1          0.7701     1.2986    0.5188    1.1430
## quartile1     0.4576     2.1855    0.3319    0.6308
## quartile2     0.9423     1.0612    0.7128    1.2457
## quartile3     0.8996     1.1115    0.6787    1.1925
## 
## Concordance= 0.701  (se = 0.016 )
## Rsquare= 0.271   (max possible= 1 )
## Likelihood ratio test= 176.4  on 10 df,   p=<2e-16
## Wald test            = 162.9  on 10 df,   p=<2e-16
## Score (logrank) test = 173.4  on 10 df,   p=<2e-16

#Result: Only Q1 is significant

contrasts(imgs$quartile) <- contr.treatment(4, base=3)
cox_qt_base3 <- coxph(Surv(survival, censor) ~ Age+Male+Immsupp+non_white+APACHE+copd+chf+quartile, data=imgs)
summary(cox_qt_base3)

## Call:
## coxph(formula = Surv(survival, censor) ~ Age + Male + Immsupp + 
##     non_white + APACHE + copd + chf + quartile, data = imgs)
## 
##   n= 559, number of events= 369 
##    (1 observation deleted due to missingness)
## 
##                 coef exp(coef)  se(coef)      z Pr(>|z|)    
## Age         0.015858  1.015984  0.003519  4.506 6.60e-06 ***
## Male1       0.054028  1.055514  0.106463  0.507 0.611817    
## Immsupp1    0.995257  2.705419  0.111847  8.898  < 2e-16 ***
## non_white1 -0.347144  0.706703  0.141602 -2.452 0.014224 *  
## APACHE      0.025299  1.025622  0.006914  3.659 0.000253 ***
## copd1       0.052707  1.054120  0.147045  0.358 0.720015    
## chf1       -0.261249  0.770089  0.201493 -1.297 0.194780    
## quartile1  -0.676079  0.508607  0.165754 -4.079 4.53e-05 ***
## quartile2   0.046316  1.047405  0.147097  0.315 0.752863    
## quartile4   0.105750  1.111544  0.143789  0.735 0.462063    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##            exp(coef) exp(-coef) lower .95 upper .95
## Age           1.0160     0.9843    1.0090    1.0230
## Male1         1.0555     0.9474    0.8567    1.3004
## Immsupp1      2.7054     0.3696    2.1729    3.3685
## non_white1    0.7067     1.4150    0.5354    0.9328
## APACHE        1.0256     0.9750    1.0118    1.0396
## copd1         1.0541     0.9487    0.7902    1.4062
## chf1          0.7701     1.2986    0.5188    1.1430
## quartile1     0.5086     1.9662    0.3675    0.7038
## quartile2     1.0474     0.9547    0.7851    1.3974
## quartile4     1.1115     0.8996    0.8386    1.4734
## 
## Concordance= 0.701  (se = 0.016 )
## Rsquare= 0.271   (max possible= 1 )
## Likelihood ratio test= 176.4  on 10 df,   p=<2e-16
## Wald test            = 162.9  on 10 df,   p=<2e-16
## Score (logrank) test = 173.4  on 10 df,   p=<2e-16

#Result: Only Q1 is significant

contrasts(imgs$quartile) <- contr.treatment(4, base=2)
cox_qt_base2 <- coxph(Surv(survival, censor) ~ Age+Male+Immsupp+non_white+APACHE+copd+chf+quartile, data=imgs)
summary(cox_qt_base2)

## Call:
## coxph(formula = Surv(survival, censor) ~ Age + Male + Immsupp + 
##     non_white + APACHE + copd + chf + quartile, data = imgs)
## 
##   n= 559, number of events= 369 
##    (1 observation deleted due to missingness)
## 
##                 coef exp(coef)  se(coef)      z Pr(>|z|)    
## Age         0.015858  1.015984  0.003519  4.506 6.60e-06 ***
## Male1       0.054028  1.055514  0.106463  0.507 0.611817    
## Immsupp1    0.995257  2.705419  0.111847  8.898  < 2e-16 ***
## non_white1 -0.347144  0.706703  0.141602 -2.452 0.014224 *  
## APACHE      0.025299  1.025622  0.006914  3.659 0.000253 ***
## copd1       0.052707  1.054120  0.147045  0.358 0.720015    
## chf1       -0.261249  0.770089  0.201493 -1.297 0.194780    
## quartile1  -0.722395  0.485588  0.163848 -4.409 1.04e-05 ***
## quartile3  -0.046316  0.954740  0.147097 -0.315 0.752863    
## quartile4   0.059434  1.061236  0.142399  0.417 0.676403    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##            exp(coef) exp(-coef) lower .95 upper .95
## Age           1.0160     0.9843    1.0090    1.0230
## Male1         1.0555     0.9474    0.8567    1.3004
## Immsupp1      2.7054     0.3696    2.1729    3.3685
## non_white1    0.7067     1.4150    0.5354    0.9328
## APACHE        1.0256     0.9750    1.0118    1.0396
## copd1         1.0541     0.9487    0.7902    1.4062
## chf1          0.7701     1.2986    0.5188    1.1430
## quartile1     0.4856     2.0594    0.3522    0.6695
## quartile3     0.9547     1.0474    0.7156    1.2738
## quartile4     1.0612     0.9423    0.8028    1.4029
## 
## Concordance= 0.701  (se = 0.016 )
## Rsquare= 0.271   (max possible= 1 )
## Likelihood ratio test= 176.4  on 10 df,   p=<2e-16
## Wald test            = 162.9  on 10 df,   p=<2e-16
## Score (logrank) test = 173.4  on 10 df,   p=<2e-16

#Result: Only Q1 is significant

#reset factor to base 1
contrasts(imgs$quartile) <- contr.treatment(4, base=1)
cox_qt_base1 <- coxph(Surv(survival, censor) ~ Age+Male+Immsupp+non_white+APACHE+copd+chf+quartile, data=imgs)
summary(cox_qt_base1)

## Call:
## coxph(formula = Surv(survival, censor) ~ Age + Male + Immsupp + 
##     non_white + APACHE + copd + chf + quartile, data = imgs)
## 
##   n= 559, number of events= 369 
##    (1 observation deleted due to missingness)
## 
##                 coef exp(coef)  se(coef)      z Pr(>|z|)    
## Age         0.015858  1.015984  0.003519  4.506 6.60e-06 ***
## Male1       0.054028  1.055514  0.106463  0.507 0.611817    
## Immsupp1    0.995257  2.705419  0.111847  8.898  < 2e-16 ***
## non_white1 -0.347144  0.706703  0.141602 -2.452 0.014224 *  
## APACHE      0.025299  1.025622  0.006914  3.659 0.000253 ***
## copd1       0.052707  1.054120  0.147045  0.358 0.720015    
## chf1       -0.261249  0.770089  0.201493 -1.297 0.194780    
## quartile2   0.722395  2.059360  0.163848  4.409 1.04e-05 ***
## quartile3   0.676079  1.966154  0.165754  4.079 4.53e-05 ***
## quartile4   0.781829  2.185466  0.163799  4.773 1.81e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##            exp(coef) exp(-coef) lower .95 upper .95
## Age           1.0160     0.9843    1.0090    1.0230
## Male1         1.0555     0.9474    0.8567    1.3004
## Immsupp1      2.7054     0.3696    2.1729    3.3685
## non_white1    0.7067     1.4150    0.5354    0.9328
## APACHE        1.0256     0.9750    1.0118    1.0396
## copd1         1.0541     0.9487    0.7902    1.4062
## chf1          0.7701     1.2986    0.5188    1.1430
## quartile2     2.0594     0.4856    1.4937    2.8392
## quartile3     1.9662     0.5086    1.4208    2.7209
## quartile4     2.1855     0.4576    1.5853    3.0128
## 
## Concordance= 0.701  (se = 0.016 )
## Rsquare= 0.271   (max possible= 1 )
## Likelihood ratio test= 176.4  on 10 df,   p=<2e-16
## Wald test            = 162.9  on 10 df,   p=<2e-16
## Score (logrank) test = 173.4  on 10 df,   p=<2e-16

plot_all <- ggsurvplot(cox_qt, data=imgs, pval=TRUE, pval.size=3,
           ggtheme = theme_minimal(), legend.title="Quartile", legend.labs=c("Q1","Q2","Q3","Q4"),
           risk.table=TRUE, tables.height = 0.2, tables.theme = theme_cleantable(),
           xlim=c(0,3650), break.time.by=365.25, xscale="d_m", ylab="",
           xlab="Months", font.legend=8, title="Entire cohort", font.title=10,
           legend=c(1.1,0.6), risk.table.y.text=FALSE,
           risk.table.fontsize=3, risk.table.title="") 
plot_all$plot <- plot_all$plot + theme(plot.title = element_text(hjust = 0.5))
plot_all

Cox regression for non-ARDS subset

#non-ARDS patients
imgs_non <- imgs %>% filter(Final_Dx != "ARDS" & Final_Dx != "Sepsis/ARDS") %>% mutate(quartile = as.factor(ntile(Total,4)))

cox_all_non <- coxph(Surv(survival, censor) ~ Age+Male+Immsupp+non_white+APACHE+copd+chf+Total, data=imgs_non) #for betas
summary(cox_all_non)

## Call:
## coxph(formula = Surv(survival, censor) ~ Age + Male + Immsupp + 
##     non_white + APACHE + copd + chf + Total, data = imgs_non)
## 
##   n= 436, number of events= 280 
##    (1 observation deleted due to missingness)
## 
##                 coef exp(coef)  se(coef)      z Pr(>|z|)    
## Age         0.019080  1.019263  0.003961  4.816 1.46e-06 ***
## Male1       0.082921  1.086456  0.123718  0.670  0.50270    
## Immsupp1    0.744200  2.104758  0.129506  5.746 9.11e-09 ***
## non_white1 -0.409955  0.663680  0.157584 -2.602  0.00928 ** 
## APACHE      0.017896  1.018057  0.008259  2.167  0.03025 *  
## copd1       0.030619  1.031093  0.160688  0.191  0.84888    
## chf1       -0.429238  0.651005  0.227067 -1.890  0.05871 .  
## Total       0.068074  1.070444  0.014178  4.801 1.58e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##            exp(coef) exp(-coef) lower .95 upper .95
## Age           1.0193     0.9811    1.0114    1.0272
## Male1         1.0865     0.9204    0.8525    1.3846
## Immsupp1      2.1048     0.4751    1.6329    2.7129
## non_white1    0.6637     1.5067    0.4873    0.9038
## APACHE        1.0181     0.9823    1.0017    1.0347
## copd1         1.0311     0.9698    0.7525    1.4128
## chf1          0.6510     1.5361    0.4172    1.0159
## Total         1.0704     0.9342    1.0411    1.1006
## 
## Concordance= 0.68  (se = 0.018 )
## Rsquare= 0.225   (max possible= 0.999 )
## Likelihood ratio test= 111.3  on 8 df,   p=<2e-16
## Wald test            = 106.5  on 8 df,   p=<2e-16
## Score (logrank) test = 111.1  on 8 df,   p=<2e-16

#std errors checked -> ok

cox_all_non_strata <- survfit(Surv(survival, censor) ~ quartile, data = imgs_non) #for KM

#std errors checked ->okay

contrasts(imgs_non$quartile) <- contr.treatment(4, base=4)
cox_non_base4 <- coxph(Surv(survival, censor) ~ Age+Male+Immsupp+non_white+APACHE+copd+chf+quartile, data=imgs_non)
summary(cox_non_base4)

## Call:
## coxph(formula = Surv(survival, censor) ~ Age + Male + Immsupp + 
##     non_white + APACHE + copd + chf + quartile, data = imgs_non)
## 
##   n= 436, number of events= 280 
##    (1 observation deleted due to missingness)
## 
##                 coef exp(coef)  se(coef)      z Pr(>|z|)    
## Age         0.017938  1.018100  0.003939  4.553 5.28e-06 ***
## Male1       0.083598  1.087192  0.123690  0.676  0.49912    
## Immsupp1    0.791019  2.205644  0.130378  6.067 1.30e-09 ***
## non_white1 -0.416692  0.659224  0.157258 -2.650  0.00806 ** 
## APACHE      0.015815  1.015940  0.008309  1.903  0.05698 .  
## copd1       0.023779  1.024064  0.163289  0.146  0.88422    
## chf1       -0.396363  0.672763  0.229551 -1.727  0.08422 .  
## quartile1  -0.918378  0.399166  0.191878 -4.786 1.70e-06 ***
## quartile2  -0.330150  0.718816  0.170259 -1.939  0.05249 .  
## quartile3  -0.079902  0.923207  0.161820 -0.494  0.62147    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##            exp(coef) exp(-coef) lower .95 upper .95
## Age           1.0181     0.9822    1.0103    1.0260
## Male1         1.0872     0.9198    0.8531    1.3855
## Immsupp1      2.2056     0.4534    1.7083    2.8478
## non_white1    0.6592     1.5169    0.4844    0.8972
## APACHE        1.0159     0.9843    0.9995    1.0326
## copd1         1.0241     0.9765    0.7436    1.4103
## chf1          0.6728     1.4864    0.4290    1.0550
## quartile1     0.3992     2.5052    0.2740    0.5814
## quartile2     0.7188     1.3912    0.5149    1.0036
## quartile3     0.9232     1.0832    0.6723    1.2678
## 
## Concordance= 0.684  (se = 0.018 )
## Rsquare= 0.238   (max possible= 0.999 )
## Likelihood ratio test= 118.4  on 10 df,   p=<2e-16
## Wald test            = 109.5  on 10 df,   p=<2e-16
## Score (logrank) test = 115.3  on 10 df,   p=<2e-16

#Result: Only Q1 is significant

contrasts(imgs_non$quartile) <- contr.treatment(4, base=3)
cox_non_base3 <- coxph(Surv(survival, censor) ~ Age+Male+Immsupp+non_white+APACHE+copd+chf+quartile, data=imgs_non)
summary(cox_non_base3)

## Call:
## coxph(formula = Surv(survival, censor) ~ Age + Male + Immsupp + 
##     non_white + APACHE + copd + chf + quartile, data = imgs_non)
## 
##   n= 436, number of events= 280 
##    (1 observation deleted due to missingness)
## 
##                 coef exp(coef)  se(coef)      z Pr(>|z|)    
## Age         0.017938  1.018100  0.003939  4.553 5.28e-06 ***
## Male1       0.083598  1.087192  0.123690  0.676  0.49912    
## Immsupp1    0.791019  2.205644  0.130378  6.067 1.30e-09 ***
## non_white1 -0.416692  0.659224  0.157258 -2.650  0.00806 ** 
## APACHE      0.015815  1.015940  0.008309  1.903  0.05698 .  
## copd1       0.023779  1.024064  0.163289  0.146  0.88422    
## chf1       -0.396363  0.672763  0.229551 -1.727  0.08422 .  
## quartile1  -0.838476  0.432369  0.186605 -4.493 7.01e-06 ***
## quartile2  -0.250249  0.778607  0.166932 -1.499  0.13385    
## quartile4   0.079902  1.083180  0.161820  0.494  0.62147    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##            exp(coef) exp(-coef) lower .95 upper .95
## Age           1.0181     0.9822    1.0103    1.0260
## Male1         1.0872     0.9198    0.8531    1.3855
## Immsupp1      2.2056     0.4534    1.7083    2.8478
## non_white1    0.6592     1.5169    0.4844    0.8972
## APACHE        1.0159     0.9843    0.9995    1.0326
## copd1         1.0241     0.9765    0.7436    1.4103
## chf1          0.6728     1.4864    0.4290    1.0550
## quartile1     0.4324     2.3128    0.2999    0.6233
## quartile2     0.7786     1.2843    0.5613    1.0800
## quartile4     1.0832     0.9232    0.7888    1.4874
## 
## Concordance= 0.684  (se = 0.018 )
## Rsquare= 0.238   (max possible= 0.999 )
## Likelihood ratio test= 118.4  on 10 df,   p=<2e-16
## Wald test            = 109.5  on 10 df,   p=<2e-16
## Score (logrank) test = 115.3  on 10 df,   p=<2e-16

#Result: Only Q1 is significant

contrasts(imgs_non$quartile) <- contr.treatment(4, base=2)
cox_non_base2 <- coxph(Surv(survival, censor) ~ Age+Male+Immsupp+non_white+APACHE+copd+chf+quartile, data=imgs_non)
summary(cox_non_base2)

## Call:
## coxph(formula = Surv(survival, censor) ~ Age + Male + Immsupp + 
##     non_white + APACHE + copd + chf + quartile, data = imgs_non)
## 
##   n= 436, number of events= 280 
##    (1 observation deleted due to missingness)
## 
##                 coef exp(coef)  se(coef)      z Pr(>|z|)    
## Age         0.017938  1.018100  0.003939  4.553 5.28e-06 ***
## Male1       0.083598  1.087192  0.123690  0.676  0.49912    
## Immsupp1    0.791019  2.205644  0.130378  6.067 1.30e-09 ***
## non_white1 -0.416692  0.659224  0.157258 -2.650  0.00806 ** 
## APACHE      0.015815  1.015940  0.008309  1.903  0.05698 .  
## copd1       0.023779  1.024064  0.163289  0.146  0.88422    
## chf1       -0.396363  0.672763  0.229551 -1.727  0.08422 .  
## quartile1  -0.588228  0.555311  0.194310 -3.027  0.00247 ** 
## quartile3   0.250249  1.284345  0.166932  1.499  0.13385    
## quartile4   0.330150  1.391177  0.170259  1.939  0.05249 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##            exp(coef) exp(-coef) lower .95 upper .95
## Age           1.0181     0.9822    1.0103    1.0260
## Male1         1.0872     0.9198    0.8531    1.3855
## Immsupp1      2.2056     0.4534    1.7083    2.8478
## non_white1    0.6592     1.5169    0.4844    0.8972
## APACHE        1.0159     0.9843    0.9995    1.0326
## copd1         1.0241     0.9765    0.7436    1.4103
## chf1          0.6728     1.4864    0.4290    1.0550
## quartile1     0.5553     1.8008    0.3794    0.8127
## quartile3     1.2843     0.7786    0.9260    1.7815
## quartile4     1.3912     0.7188    0.9965    1.9423
## 
## Concordance= 0.684  (se = 0.018 )
## Rsquare= 0.238   (max possible= 0.999 )
## Likelihood ratio test= 118.4  on 10 df,   p=<2e-16
## Wald test            = 109.5  on 10 df,   p=<2e-16
## Score (logrank) test = 115.3  on 10 df,   p=<2e-16

#Result: Only Q1 is significant

#reset factor to base 1
contrasts(imgs_non$quartile) <- contr.treatment(4, base=1)
cox_non_base1 <- coxph(Surv(survival, censor) ~ Age+Male+Immsupp+non_white+APACHE+copd+chf+quartile, data=imgs_non)
summary(cox_non_base1)

## Call:
## coxph(formula = Surv(survival, censor) ~ Age + Male + Immsupp + 
##     non_white + APACHE + copd + chf + quartile, data = imgs_non)
## 
##   n= 436, number of events= 280 
##    (1 observation deleted due to missingness)
## 
##                 coef exp(coef)  se(coef)      z Pr(>|z|)    
## Age         0.017938  1.018100  0.003939  4.553 5.28e-06 ***
## Male1       0.083598  1.087192  0.123690  0.676  0.49912    
## Immsupp1    0.791019  2.205644  0.130378  6.067 1.30e-09 ***
## non_white1 -0.416692  0.659224  0.157258 -2.650  0.00806 ** 
## APACHE      0.015815  1.015940  0.008309  1.903  0.05698 .  
## copd1       0.023779  1.024064  0.163289  0.146  0.88422    
## chf1       -0.396363  0.672763  0.229551 -1.727  0.08422 .  
## quartile2   0.588228  1.800794  0.194310  3.027  0.00247 ** 
## quartile3   0.838476  2.312840  0.186605  4.493 7.01e-06 ***
## quartile4   0.918378  2.505223  0.191878  4.786 1.70e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##            exp(coef) exp(-coef) lower .95 upper .95
## Age           1.0181     0.9822    1.0103    1.0260
## Male1         1.0872     0.9198    0.8531    1.3855
## Immsupp1      2.2056     0.4534    1.7083    2.8478
## non_white1    0.6592     1.5169    0.4844    0.8972
## APACHE        1.0159     0.9843    0.9995    1.0326
## copd1         1.0241     0.9765    0.7436    1.4103
## chf1          0.6728     1.4864    0.4290    1.0550
## quartile2     1.8008     0.5553    1.2305    2.6355
## quartile3     2.3128     0.4324    1.6044    3.3341
## quartile4     2.5052     0.3992    1.7200    3.6490
## 
## Concordance= 0.684  (se = 0.018 )
## Rsquare= 0.238   (max possible= 0.999 )
## Likelihood ratio test= 118.4  on 10 df,   p=<2e-16
## Wald test            = 109.5  on 10 df,   p=<2e-16
## Score (logrank) test = 115.3  on 10 df,   p=<2e-16

plot_non <- ggsurvplot(cox_all_non_strata, data=imgs_non, pval=TRUE, pval.size=3,
           ggtheme = theme_minimal(), legend.title="Quartile", 
           risk.table=TRUE, tables.height = 0.2, tables.theme = theme_cleantable(),
           xlim=c(0,3650), break.time.by=365.25, xscale="d_m", ylab="",
           xlab="Months", font.legend=8, title="Non-ARDS patients", font.title=10,
           legend="none", risk.table.y.text=FALSE, 
           risk.table.fontsize=3, risk.table.title="") 
plot_non$plot <- plot_non$plot + theme(axis.text.y = element_blank(), plot.title = element_text(hjust = 0.5))
plot_non

KM curves side by side

KM_curves <- list()
KM_curves[[1]] <- plot_all
KM_curves[[2]] <- plot_non
arrange_ggsurvplots(KM_curves, print=TRUE, ncol=2, nrow=1)

#std errors checked ->okay
hosp_mort <- glm(hosp_death ~ quartile+Male+Age+Immsupp+APACHE+non_white+copd+chf, family=binomial(link='logit'), data=imgs)
summary(hosp_mort)

## 
## Call:
## glm(formula = hosp_death ~ quartile + Male + Age + Immsupp + 
##     APACHE + non_white + copd + chf, family = binomial(link = "logit"), 
##     data = imgs)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5878  -0.7771  -0.4981   0.8897   2.5633  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -4.308019   0.619638  -6.952 3.59e-12 ***
## quartile2    0.977419   0.350667   2.787 0.005315 ** 
## quartile3    1.273048   0.344812   3.692 0.000222 ***
## quartile4    1.313196   0.342296   3.836 0.000125 ***
## Male1       -0.007186   0.211581  -0.034 0.972908    
## Age          0.007612   0.006892   1.104 0.269390    
## Immsupp1     1.284147   0.215892   5.948 2.71e-09 ***
## APACHE       0.058885   0.013126   4.486 7.25e-06 ***
## non_white1  -0.533414   0.294295  -1.813 0.069907 .  
## copd1       -0.132786   0.302055  -0.440 0.660220    
## chf1        -0.728473   0.417400  -1.745 0.080939 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 667.55  on 558  degrees of freedom
## Residual deviance: 556.32  on 548  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: 578.32
## 
## Number of Fisher Scoring iterations: 5

cbind(OR = exp(coef(hosp_mort)),
      exp(confint(hosp_mort)),
      p = coef(summary(hosp_mort))[,4])

##                     OR       2.5 %     97.5 %            p
## (Intercept) 0.01346018 0.003827569 0.04361803 3.589179e-12
## quartile2   2.65758874 1.357107104 5.40098142 5.314693e-03
## quartile3   3.57172343 1.848338033 7.18684749 2.224931e-04
## quartile4   3.71803703 1.934507631 7.44744932 1.248334e-04
## Male1       0.99284009 0.655723601 1.50471217 9.729076e-01
## Age         1.00764099 0.994182653 1.02146055 2.693897e-01
## Immsupp1    3.61158417 2.374076826 5.54060270 2.712762e-09
## APACHE      1.06065325 1.034088295 1.08880700 7.247657e-06
## non_white1  0.58659886 0.322492002 1.02783849 6.990708e-02
## copd1       0.87565251 0.477472759 1.56698187 6.602203e-01
## chf1        0.48264533 0.202196883 1.05496607 8.093892e-02

contrasts(imgs$quartile) <- contr.treatment(4, base=2)
inhosp_quart_base2 <- glm(hosp_death ~ quartile+Male+Age+Immsupp+APACHE+non_white+copd+chf, family=binomial(link='logit'), data=imgs)
summary(inhosp_quart_base2)

## 
## Call:
## glm(formula = hosp_death ~ quartile + Male + Age + Immsupp + 
##     APACHE + non_white + copd + chf, family = binomial(link = "logit"), 
##     data = imgs)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5878  -0.7771  -0.4981   0.8897   2.5633  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.330600   0.618051  -5.389 7.09e-08 ***
## quartile1   -0.977419   0.350667  -2.787  0.00531 ** 
## quartile3    0.295629   0.283931   1.041  0.29778    
## quartile4    0.335777   0.279645   1.201  0.22986    
## Male1       -0.007186   0.211581  -0.034  0.97291    
## Age          0.007612   0.006892   1.104  0.26939    
## Immsupp1     1.284147   0.215892   5.948 2.71e-09 ***
## APACHE       0.058885   0.013126   4.486 7.25e-06 ***
## non_white1  -0.533414   0.294295  -1.813  0.06991 .  
## copd1       -0.132786   0.302055  -0.440  0.66022    
## chf1        -0.728473   0.417400  -1.745  0.08094 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 667.55  on 558  degrees of freedom
## Residual deviance: 556.32  on 548  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: 578.32
## 
## Number of Fisher Scoring iterations: 5

#result: Q1is significant, 3 and 4 are not

cbind(OR = exp(coef(inhosp_quart_base2)),
      exp(confint(inhosp_quart_base2)),
      p = coef(summary(inhosp_quart_base2))[,4])

##                     OR      2.5 %    97.5 %            p
## (Intercept) 0.03577163 0.01031823 0.1168662 7.089952e-08
## quartile1   0.37628094 0.18515153 0.7368615 5.314693e-03
## quartile3   1.34397147 0.77105057 2.3519011 2.977824e-01
## quartile4   1.39902649 0.80976693 2.4284928 2.298585e-01
## Male1       0.99284009 0.65572360 1.5047122 9.729076e-01
## Age         1.00764099 0.99418265 1.0214605 2.693897e-01
## Immsupp1    3.61158417 2.37407683 5.5406027 2.712762e-09
## APACHE      1.06065325 1.03408829 1.0888070 7.247657e-06
## non_white1  0.58659886 0.32249200 1.0278385 6.990708e-02
## copd1       0.87565251 0.47747276 1.5669819 6.602203e-01
## chf1        0.48264533 0.20219688 1.0549661 8.093892e-02

contrasts(imgs$quartile) <- contr.treatment(4, base=3)
inhosp_quart_base3 <- glm(hosp_death ~ quartile+Male+Age+Immsupp+APACHE+non_white+copd+chf, family=binomial(link='logit'), data=imgs)
summary(inhosp_quart_base3)

## 
## Call:
## glm(formula = hosp_death ~ quartile + Male + Age + Immsupp + 
##     APACHE + non_white + copd + chf, family = binomial(link = "logit"), 
##     data = imgs)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5878  -0.7771  -0.4981   0.8897   2.5633  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.034971   0.599423  -5.063 4.12e-07 ***
## quartile1   -1.273048   0.344812  -3.692 0.000222 ***
## quartile2   -0.295629   0.283931  -1.041 0.297782    
## quartile4    0.040148   0.271342   0.148 0.882375    
## Male1       -0.007186   0.211581  -0.034 0.972908    
## Age          0.007612   0.006892   1.104 0.269390    
## Immsupp1     1.284147   0.215892   5.948 2.71e-09 ***
## APACHE       0.058885   0.013126   4.486 7.25e-06 ***
## non_white1  -0.533414   0.294295  -1.813 0.069907 .  
## copd1       -0.132786   0.302055  -0.440 0.660220    
## chf1        -0.728473   0.417400  -1.745 0.080939 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 667.55  on 558  degrees of freedom
## Residual deviance: 556.32  on 548  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: 578.32
## 
## Number of Fisher Scoring iterations: 5

#result: Q1 is significant, 2 and 4 are not

cbind(OR = exp(coef(inhosp_quart_base3)),
      exp(confint(inhosp_quart_base3)),
      p = coef(summary(inhosp_quart_base3))[,4])

##                     OR      2.5 %    97.5 %            p
## (Intercept) 0.04807605 0.01441297 0.1517703 4.123787e-07
## quartile1   0.27997689 0.13914307 0.5410266 2.224931e-04
## quartile2   0.74406342 0.42518794 1.2969318 2.977824e-01
## quartile4   1.04096443 0.61117357 1.7737158 8.823749e-01
## Male1       0.99284009 0.65572360 1.5047122 9.729076e-01
## Age         1.00764099 0.99418265 1.0214605 2.693897e-01
## Immsupp1    3.61158417 2.37407683 5.5406027 2.712762e-09
## APACHE      1.06065325 1.03408829 1.0888070 7.247657e-06
## non_white1  0.58659886 0.32249200 1.0278385 6.990708e-02
## copd1       0.87565251 0.47747276 1.5669819 6.602203e-01
## chf1        0.48264533 0.20219688 1.0549661 8.093892e-02

contrasts(imgs$quartile) <- contr.treatment(4, base=4)
inhosp_quart_base4 <- glm(hosp_death ~ quartile+Male+Age+Immsupp+APACHE+non_white+copd+chf, family=binomial(link='logit'), data=imgs)
summary(inhosp_quart_base4)

## 
## Call:
## glm(formula = hosp_death ~ quartile + Male + Age + Immsupp + 
##     APACHE + non_white + copd + chf, family = binomial(link = "logit"), 
##     data = imgs)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5878  -0.7771  -0.4981   0.8897   2.5633  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.994823   0.608013  -4.926 8.41e-07 ***
## quartile1   -1.313196   0.342296  -3.836 0.000125 ***
## quartile2   -0.335777   0.279645  -1.201 0.229858    
## quartile3   -0.040148   0.271342  -0.148 0.882375    
## Male1       -0.007186   0.211581  -0.034 0.972908    
## Age          0.007612   0.006892   1.104 0.269390    
## Immsupp1     1.284147   0.215892   5.948 2.71e-09 ***
## APACHE       0.058885   0.013126   4.486 7.25e-06 ***
## non_white1  -0.533414   0.294295  -1.813 0.069907 .  
## copd1       -0.132786   0.302055  -0.440 0.660220    
## chf1        -0.728473   0.417400  -1.745 0.080939 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 667.55  on 558  degrees of freedom
## Residual deviance: 556.32  on 548  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: 578.32
## 
## Number of Fisher Scoring iterations: 5

#result: Q1 is significant, 2 and 3 are not

cbind(OR = exp(coef(inhosp_quart_base4)),
      exp(confint(inhosp_quart_base4)),
      p = coef(summary(inhosp_quart_base4))[,4])

##                     OR      2.5 %    97.5 %            p
## (Intercept) 0.05004546 0.01474234 0.1605245 8.410630e-07
## quartile1   0.26895913 0.13427416 0.5169274 1.248334e-04
## quartile2   0.71478275 0.41177804 1.2349232 2.298585e-01
## quartile3   0.96064762 0.56378818 1.6361964 8.823749e-01
## Male1       0.99284009 0.65572360 1.5047122 9.729076e-01
## Age         1.00764099 0.99418265 1.0214605 2.693897e-01
## Immsupp1    3.61158417 2.37407683 5.5406027 2.712762e-09
## APACHE      1.06065325 1.03408829 1.0888070 7.247657e-06
## non_white1  0.58659886 0.32249200 1.0278385 6.990708e-02
## copd1       0.87565251 0.47747276 1.5669819 6.602203e-01
## chf1        0.48264533 0.20219688 1.0549661 8.093892e-02

#reset factor to base 1
contrasts(imgs$quartile) <- contr.treatment(4, base=1)

#std errors checked ->okay

hosp_mort_non_r1 <- glm(hosp_death ~ quartile+Male+Age+Immsupp+APACHE+non_white+copd+chf, family=binomial(link='logit'), data=imgs_non)
summary(hosp_mort_non_r1)

## 
## Call:
## glm(formula = hosp_death ~ quartile + Male + Age + Immsupp + 
##     APACHE + non_white + copd + chf, family = binomial(link = "logit"), 
##     data = imgs_non)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.4147  -0.6999  -0.4993  -0.2650   2.5139  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.656926   0.695942  -5.255 1.48e-07 ***
## quartile2    0.422174   0.414174   1.019 0.308052    
## quartile3    1.041093   0.382233   2.724 0.006455 ** 
## quartile4    0.863732   0.393273   2.196 0.028073 *  
## Male1       -0.114983   0.253055  -0.454 0.649556    
## Age          0.007204   0.008351   0.863 0.388323    
## Immsupp1     0.918398   0.259996   3.532 0.000412 ***
## APACHE       0.047948   0.016177   2.964 0.003037 ** 
## non_white1  -0.843357   0.373297  -2.259 0.023870 *  
## copd1       -0.123341   0.346305  -0.356 0.721718    
## chf1        -0.862727   0.525182  -1.643 0.100441    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 451.96  on 435  degrees of freedom
## Residual deviance: 398.98  on 425  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: 420.98
## 
## Number of Fisher Scoring iterations: 5

cbind(OR = exp(coef(hosp_mort_non_r1)),
      exp(confint(hosp_mort_non_r1)),
      p = coef(summary(hosp_mort_non_r1))[,4])

##                     OR       2.5 %     97.5 %            p
## (Intercept) 0.02581175 0.006256306 0.09638051 1.483154e-07
## quartile2   1.52527456 0.682180953 3.49756666 3.080525e-01
## quartile3   2.83231015 1.363691206 6.15708813 6.455294e-03
## quartile4   2.37199718 1.113636707 5.25176753 2.807291e-02
## Male1       0.89138126 0.542383131 1.46589497 6.495556e-01
## Age         1.00723006 0.991001173 1.02407737 3.883234e-01
## Immsupp1    2.50527300 1.507060048 4.18531889 4.118774e-04
## APACHE      1.04911626 1.016565867 1.08333341 3.036661e-03
## non_white1  0.43026373 0.196597525 0.86158480 2.387012e-02
## copd1       0.88396182 0.437157094 1.71182926 7.217176e-01
## chf1        0.42200962 0.134620177 1.09600071 1.004410e-01

#ref level =2
contrasts(imgs_non$quartile) <- contr.treatment(4, base = 2)
hosp_mort_non_r2 <- glm(hosp_death ~ quartile+Male+Age+Immsupp+APACHE+non_white+copd+chf, family=binomial(link='logit'), data=imgs_non)
summary(hosp_mort_non_r2)

## 
## Call:
## glm(formula = hosp_death ~ quartile + Male + Age + Immsupp + 
##     APACHE + non_white + copd + chf, family = binomial(link = "logit"), 
##     data = imgs_non)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.4147  -0.6999  -0.4993  -0.2650   2.5139  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.234751   0.738440  -4.381 1.18e-05 ***
## quartile1   -0.422174   0.414174  -1.019 0.308052    
## quartile3    0.618918   0.349294   1.772 0.076409 .  
## quartile4    0.441558   0.358730   1.231 0.218363    
## Male1       -0.114983   0.253055  -0.454 0.649556    
## Age          0.007204   0.008351   0.863 0.388323    
## Immsupp1     0.918398   0.259996   3.532 0.000412 ***
## APACHE       0.047948   0.016177   2.964 0.003037 ** 
## non_white1  -0.843357   0.373297  -2.259 0.023870 *  
## copd1       -0.123341   0.346305  -0.356 0.721718    
## chf1        -0.862727   0.525182  -1.643 0.100441    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 451.96  on 435  degrees of freedom
## Residual deviance: 398.98  on 425  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: 420.98
## 
## Number of Fisher Scoring iterations: 5

cbind(OR = exp(coef(hosp_mort_non_r2)),
      exp(confint(hosp_mort_non_r2)),
      p = coef(summary(hosp_mort_non_r2))[,4])

##                    OR       2.5 %    97.5 %            p
## (Intercept) 0.0393700 0.008820534 0.1606021 1.183977e-05
## quartile1   0.6556197 0.285913064 1.4658867 3.080525e-01
## quartile3   1.8569182 0.942948850 3.7287019 7.640892e-02
## quartile4   1.5551280 0.773205563 3.1735108 2.183629e-01
## Male1       0.8913813 0.542383131 1.4658950 6.495556e-01
## Age         1.0072301 0.991001173 1.0240774 3.883234e-01
## Immsupp1    2.5052730 1.507060048 4.1853189 4.118774e-04
## APACHE      1.0491163 1.016565867 1.0833334 3.036661e-03
## non_white1  0.4302637 0.196597525 0.8615848 2.387012e-02
## copd1       0.8839618 0.437157094 1.7118293 7.217176e-01
## chf1        0.4220096 0.134620177 1.0960007 1.004410e-01

#ref level = 3
contrasts(imgs_non$quartile) <- contr.treatment(4, base = 3)
hosp_mort_non_r3 <- glm(hosp_death ~ quartile+Male+Age+Immsupp+APACHE+non_white+copd+chf, family=binomial(link='logit'), data=imgs_non)
summary(hosp_mort_non_r3)

## 
## Call:
## glm(formula = hosp_death ~ quartile + Male + Age + Immsupp + 
##     APACHE + non_white + copd + chf, family = binomial(link = "logit"), 
##     data = imgs_non)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.4147  -0.6999  -0.4993  -0.2650   2.5139  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.615833   0.697808  -3.749 0.000178 ***
## quartile1   -1.041093   0.382233  -2.724 0.006455 ** 
## quartile2   -0.618918   0.349294  -1.772 0.076409 .  
## quartile4   -0.177360   0.323987  -0.547 0.584083    
## Male1       -0.114983   0.253055  -0.454 0.649556    
## Age          0.007204   0.008351   0.863 0.388323    
## Immsupp1     0.918398   0.259996   3.532 0.000412 ***
## APACHE       0.047948   0.016177   2.964 0.003037 ** 
## non_white1  -0.843357   0.373297  -2.259 0.023870 *  
## copd1       -0.123341   0.346305  -0.356 0.721718    
## chf1        -0.862727   0.525182  -1.643 0.100441    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 451.96  on 435  degrees of freedom
## Residual deviance: 398.98  on 425  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: 420.98
## 
## Number of Fisher Scoring iterations: 5

cbind(OR = exp(coef(hosp_mort_non_r3)),
      exp(confint(hosp_mort_non_r3)),
      p = coef(summary(hosp_mort_non_r3))[,4])

##                     OR      2.5 %    97.5 %            p
## (Intercept) 0.07310687 0.01786178 0.2773064 0.0001777939
## quartile1   0.35306868 0.16241444 0.7333038 0.0064552941
## quartile2   0.53852667 0.26818985 1.0605029 0.0764089248
## quartile4   0.83747791 0.44198820 1.5796181 0.5840826873
## Male1       0.89138126 0.54238313 1.4658950 0.6495556183
## Age         1.00723006 0.99100117 1.0240774 0.3883234008
## Immsupp1    2.50527300 1.50706005 4.1853189 0.0004118774
## APACHE      1.04911626 1.01656587 1.0833334 0.0030366608
## non_white1  0.43026373 0.19659752 0.8615848 0.0238701170
## copd1       0.88396182 0.43715709 1.7118293 0.7217175667
## chf1        0.42200962 0.13462018 1.0960007 0.1004410074

#ref level = 4
contrasts(imgs_non$quartile) <- contr.treatment(4, base = 4)
hosp_mort_non_r4 <- glm(hosp_death ~ quartile+Male+Age+Immsupp+APACHE+non_white+copd+chf, family=binomial(link='logit'), data=imgs_non)
summary(hosp_mort_non_r4)

## 
## Call:
## glm(formula = hosp_death ~ quartile + Male + Age + Immsupp + 
##     APACHE + non_white + copd + chf, family = binomial(link = "logit"), 
##     data = imgs_non)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.4147  -0.6999  -0.4993  -0.2650   2.5139  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.793193   0.715721  -3.903 9.52e-05 ***
## quartile1   -0.863732   0.393273  -2.196 0.028073 *  
## quartile2   -0.441558   0.358730  -1.231 0.218363    
## quartile3    0.177360   0.323987   0.547 0.584083    
## Male1       -0.114983   0.253055  -0.454 0.649556    
## Age          0.007204   0.008351   0.863 0.388323    
## Immsupp1     0.918398   0.259996   3.532 0.000412 ***
## APACHE       0.047948   0.016177   2.964 0.003037 ** 
## non_white1  -0.843357   0.373297  -2.259 0.023870 *  
## copd1       -0.123341   0.346305  -0.356 0.721718    
## chf1        -0.862727   0.525182  -1.643 0.100441    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 451.96  on 435  degrees of freedom
## Residual deviance: 398.98  on 425  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: 420.98
## 
## Number of Fisher Scoring iterations: 5

cbind(OR = exp(coef(hosp_mort_non_r4)),
      exp(confint(hosp_mort_non_r4)),
      p = coef(summary(hosp_mort_non_r4))[,4])

##                     OR      2.5 %    97.5 %            p
## (Intercept) 0.06122539 0.01437719 0.2395434 9.515319e-05
## quartile1   0.42158566 0.19041208 0.8979589 2.807291e-02
## quartile2   0.64303388 0.31510843 1.2933171 2.183629e-01
## quartile3   1.19406135 0.63306442 2.2625038 5.840827e-01
## Male1       0.89138126 0.54238313 1.4658950 6.495556e-01
## Age         1.00723006 0.99100117 1.0240774 3.883234e-01
## Immsupp1    2.50527300 1.50706005 4.1853189 4.118774e-04
## APACHE      1.04911626 1.01656587 1.0833334 3.036661e-03
## non_white1  0.43026373 0.19659752 0.8615848 2.387012e-02
## copd1       0.88396182 0.43715709 1.7118293 7.217176e-01
## chf1        0.42200962 0.13462018 1.0960007 1.004410e-01

#reset ref level to 1
contrasts(imgs_non$quartile) <- contr.treatment(4, base = 1)

ICU Length of stay

contrasts(imgs$quartile) <-contr.treatment(4,base = 1)

#Negative binomial model
#std errors checked -> okay

ICUdays_nb <- glm.nb(ICU_days ~ Total+APACHE+Age+non_white+Male+Immsupp+copd+chf, data=imgs)
summary(ICUdays_nb)

## 
## Call:
## glm.nb(formula = ICU_days ~ Total + APACHE + Age + non_white + 
##     Male + Immsupp + copd + chf, data = imgs, init.theta = 3.02729853, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.2276  -0.8067  -0.3650   0.2561   3.6493  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.493862   0.148292  10.074  < 2e-16 ***
## Total        0.056503   0.006971   8.106 5.25e-16 ***
## APACHE       0.022434   0.003631   6.178 6.48e-10 ***
## Age         -0.009249   0.001905  -4.856 1.20e-06 ***
## non_white1  -0.242192   0.075780  -3.196  0.00139 ** 
## Male1        0.100367   0.059707   1.681  0.09277 .  
## Immsupp1    -0.074787   0.063722  -1.174  0.24054    
## copd1       -0.014389   0.084868  -0.170  0.86536    
## chf1        -0.149653   0.110395  -1.356  0.17522    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(3.0273) family taken to be 1)
## 
##     Null deviance: 690.58  on 558  degrees of freedom
## Residual deviance: 533.76  on 550  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: 3117.2
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  3.027 
##           Std. Err.:  0.247 
## 
##  2 x log-likelihood:  -3097.171

cbind(IRR = exp(coef(ICUdays_nb)), exp(confint(ICUdays_nb)))[2,]

##      IRR    2.5 %   97.5 % 
## 1.058129 1.043567 1.072946

#standard errors checked -> okay
#Base = 1
contrasts(imgs$quartile) <-contr.treatment(4,base = 1)
ICUdays_nb_q1 <- glm.nb(ICU_days ~ quartile+APACHE+Age+non_white+Male+Immsupp+copd+chf, data=imgs)
summary(ICUdays_nb_q1)

## 
## Call:
## glm.nb(formula = ICU_days ~ quartile + APACHE + Age + non_white + 
##     Male + Immsupp + copd + chf, data = imgs, init.theta = 3.025792556, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3082  -0.8335  -0.3579   0.2777   3.8017  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.516790   0.148598  10.207  < 2e-16 ***
## quartile2    0.334510   0.089146   3.752 0.000175 ***
## quartile3    0.582179   0.087980   6.617 3.66e-11 ***
## quartile4    0.654711   0.087890   7.449 9.39e-14 ***
## APACHE       0.022001   0.003651   6.026 1.68e-09 ***
## Age         -0.009538   0.001923  -4.960 7.03e-07 ***
## non_white1  -0.245980   0.075787  -3.246 0.001172 ** 
## Male1        0.077991   0.059951   1.301 0.193288    
## Immsupp1    -0.076519   0.063748  -1.200 0.230009    
## copd1       -0.030197   0.085383  -0.354 0.723593    
## chf1        -0.139401   0.110515  -1.261 0.207174    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(3.0258) family taken to be 1)
## 
##     Null deviance: 690.34  on 558  degrees of freedom
## Residual deviance: 533.19  on 548  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: 3120.8
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  3.026 
##           Std. Err.:  0.247 
## 
##  2 x log-likelihood:  -3096.780

cbind(IRR = exp(coef(ICUdays_nb_q1)),
      exp(confint(ICUdays_nb_q1)),
      p = coef(summary(ICUdays_nb_q1))[,4])

##                   IRR     2.5 %    97.5 %            p
## (Intercept) 4.5575713 3.3769015 6.1566319 1.837967e-24
## quartile2   1.3972550 1.1729172 1.6648426 1.751569e-04
## quartile3   1.7899342 1.5048166 2.1296709 3.660726e-11
## quartile4   1.9245869 1.6188460 2.2888263 9.389953e-14
## APACHE      1.0222452 1.0146871 1.0298937 1.676954e-09
## Age         0.9905074 0.9866450 0.9943669 7.032483e-07
## non_white1  0.7819378 0.6742871 0.9077580 1.171743e-03
## Male1       1.0811128 0.9606191 1.2165457 1.932881e-01
## Immsupp1    0.9263357 0.8178175 1.0494715 2.300089e-01
## copd1       0.9702548 0.8213731 1.1482057 7.235930e-01
## chf1        0.8698788 0.7014168 1.0830270 2.071737e-01

#Base = 2
contrasts(imgs$quartile) <-contr.treatment(4,base = 2)
ICUdays_nb_q2 <- glm.nb(ICU_days ~ quartile+APACHE+Age+non_white+Male+Immsupp+copd+chf, data=imgs)
summary(ICUdays_nb_q2)

## 
## Call:
## glm.nb(formula = ICU_days ~ quartile + APACHE + Age + non_white + 
##     Male + Immsupp + copd + chf, data = imgs, init.theta = 3.025792556, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3082  -0.8335  -0.3579   0.2777   3.8017  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.851299   0.161818  11.441  < 2e-16 ***
## quartile1   -0.334510   0.089146  -3.752 0.000175 ***
## quartile3    0.247669   0.083791   2.956 0.003119 ** 
## quartile4    0.320202   0.083407   3.839 0.000124 ***
## APACHE       0.022001   0.003651   6.026 1.68e-09 ***
## Age         -0.009538   0.001923  -4.960 7.03e-07 ***
## non_white1  -0.245980   0.075787  -3.246 0.001172 ** 
## Male1        0.077991   0.059951   1.301 0.193288    
## Immsupp1    -0.076519   0.063748  -1.200 0.230009    
## copd1       -0.030197   0.085383  -0.354 0.723593    
## chf1        -0.139401   0.110515  -1.261 0.207174    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(3.0258) family taken to be 1)
## 
##     Null deviance: 690.34  on 558  degrees of freedom
## Residual deviance: 533.19  on 548  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: 3120.8
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  3.026 
##           Std. Err.:  0.247 
## 
##  2 x log-likelihood:  -3096.780

cbind(IRR = exp(coef(ICUdays_nb_q2)),
      exp(confint(ICUdays_nb_q2)),
      p = coef(summary(ICUdays_nb_q2))[,4])

##                   IRR     2.5 %    97.5 %            p
## (Intercept) 6.3680892 4.5987356 8.8315486 2.620331e-30
## quartile1   0.7156890 0.6006574 0.8525751 1.751569e-04
## quartile3   1.2810362 1.0871930 1.5095506 3.118715e-03
## quartile4   1.3774056 1.1684231 1.6239800 1.235251e-04
## APACHE      1.0222452 1.0146871 1.0298937 1.676954e-09
## Age         0.9905074 0.9866450 0.9943669 7.032483e-07
## non_white1  0.7819378 0.6742871 0.9077580 1.171743e-03
## Male1       1.0811128 0.9606191 1.2165457 1.932881e-01
## Immsupp1    0.9263357 0.8178175 1.0494715 2.300089e-01
## copd1       0.9702548 0.8213731 1.1482057 7.235930e-01
## chf1        0.8698788 0.7014168 1.0830270 2.071737e-01

#Base = 3
contrasts(imgs$quartile) <-contr.treatment(4,base = 3)
ICUdays_nb_q3 <- glm.nb(ICU_days ~ quartile+APACHE+Age+non_white+Male+Immsupp+copd+chf, data=imgs)
summary(ICUdays_nb_q3)

## 
## Call:
## glm.nb(formula = ICU_days ~ quartile + APACHE + Age + non_white + 
##     Male + Immsupp + copd + chf, data = imgs, init.theta = 3.025792556, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3082  -0.8335  -0.3579   0.2777   3.8017  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.098969   0.157668  13.313  < 2e-16 ***
## quartile1   -0.582179   0.087980  -6.617 3.66e-11 ***
## quartile2   -0.247669   0.083791  -2.956  0.00312 ** 
## quartile4    0.072532   0.081171   0.894  0.37155    
## APACHE       0.022001   0.003651   6.026 1.68e-09 ***
## Age         -0.009538   0.001923  -4.960 7.03e-07 ***
## non_white1  -0.245980   0.075787  -3.246  0.00117 ** 
## Male1        0.077991   0.059951   1.301  0.19329    
## Immsupp1    -0.076519   0.063748  -1.200  0.23001    
## copd1       -0.030197   0.085383  -0.354  0.72359    
## chf1        -0.139401   0.110515  -1.261  0.20717    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(3.0258) family taken to be 1)
## 
##     Null deviance: 690.34  on 558  degrees of freedom
## Residual deviance: 533.19  on 548  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: 3120.8
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  3.026 
##           Std. Err.:  0.247 
## 
##  2 x log-likelihood:  -3096.780

cbind(IRR = exp(coef(ICUdays_nb_q3)),
      exp(confint(ICUdays_nb_q3)),
      p = coef(summary(ICUdays_nb_q3))[,4])

##                   IRR     2.5 %     97.5 %            p
## (Intercept) 8.1577527 5.9060823 11.2873406 1.955984e-40
## quartile1   0.5586798 0.4695561  0.6645328 3.660726e-11
## quartile2   0.7806181 0.6624488  0.9197999 3.118715e-03
## quartile4   1.0752277 0.9167725  1.2611534 3.715475e-01
## APACHE      1.0222452 1.0146871  1.0298937 1.676954e-09
## Age         0.9905074 0.9866450  0.9943669 7.032483e-07
## non_white1  0.7819378 0.6742871  0.9077580 1.171743e-03
## Male1       1.0811128 0.9606191  1.2165457 1.932881e-01
## Immsupp1    0.9263357 0.8178175  1.0494715 2.300089e-01
## copd1       0.9702548 0.8213731  1.1482057 7.235930e-01
## chf1        0.8698788 0.7014168  1.0830270 2.071737e-01

#Base = 4
contrasts(imgs$quartile) <-contr.treatment(4,base = 4)
ICUdays_nb_q4 <- glm.nb(ICU_days ~ quartile+APACHE+Age+non_white+Male+Immsupp+copd+chf, data=imgs)
summary(ICUdays_nb_q4)

## 
## Call:
## glm.nb(formula = ICU_days ~ quartile + APACHE + Age + non_white + 
##     Male + Immsupp + copd + chf, data = imgs, init.theta = 3.025792556, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3082  -0.8335  -0.3579   0.2777   3.8017  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.171501   0.159360  13.626  < 2e-16 ***
## quartile1   -0.654711   0.087890  -7.449 9.39e-14 ***
## quartile2   -0.320202   0.083407  -3.839 0.000124 ***
## quartile3   -0.072532   0.081171  -0.894 0.371548    
## APACHE       0.022001   0.003651   6.026 1.68e-09 ***
## Age         -0.009538   0.001923  -4.960 7.03e-07 ***
## non_white1  -0.245980   0.075787  -3.246 0.001172 ** 
## Male1        0.077991   0.059951   1.301 0.193288    
## Immsupp1    -0.076519   0.063748  -1.200 0.230009    
## copd1       -0.030197   0.085383  -0.354 0.723593    
## chf1        -0.139401   0.110515  -1.261 0.207174    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(3.0258) family taken to be 1)
## 
##     Null deviance: 690.34  on 558  degrees of freedom
## Residual deviance: 533.19  on 548  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: 3120.8
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  3.026 
##           Std. Err.:  0.247 
## 
##  2 x log-likelihood:  -3096.780

cbind(IRR = exp(coef(ICUdays_nb_q4)),
      exp(confint(ICUdays_nb_q4)),
      p = coef(summary(ICUdays_nb_q4))[,4])

##                   IRR     2.5 %     97.5 %            p
## (Intercept) 8.7714418 6.3617496 12.1155580 2.790730e-42
## quartile1   0.5195920 0.4369051  0.6177240 9.389953e-14
## quartile2   0.7260026 0.6157711  0.8558543 1.235251e-04
## quartile3   0.9300355 0.7929249  1.0907832 3.715475e-01
## APACHE      1.0222452 1.0146871  1.0298937 1.676954e-09
## Age         0.9905074 0.9866450  0.9943669 7.032483e-07
## non_white1  0.7819378 0.6742871  0.9077580 1.171743e-03
## Male1       1.0811128 0.9606191  1.2165457 1.932881e-01
## Immsupp1    0.9263357 0.8178175  1.0494715 2.300089e-01
## copd1       0.9702548 0.8213731  1.1482057 7.235930e-01
## chf1        0.8698788 0.7014168  1.0830270 2.071737e-01

#Reset base level
contrasts(imgs$quartile) <-contr.treatment(4,base = 1)

Duration of mechanical ventilation

#make the subset of data: only intubated patients (n = 286)
intub <- imgs %>% filter(days_MV != 0)

#standard errors checked ->okay
MVdays_nb2 <- glm.nb(days_MV ~ Total+APACHE+Age+non_white+Male+Immsupp+copd+chf, data=intub)
summary(MVdays_nb2)

## 
## Call:
## glm.nb(formula = days_MV ~ Total + APACHE + Age + non_white + 
##     Male + Immsupp + copd + chf, data = intub, init.theta = 1.805548498, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7943  -0.9164  -0.4188   0.2213   4.3968  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.139983   0.260729   8.208 2.26e-16 ***
## Total        0.044079   0.012079   3.649 0.000263 ***
## APACHE       0.003229   0.006372   0.507 0.612389    
## Age         -0.008990   0.003171  -2.835 0.004576 ** 
## non_white1   0.016122   0.137941   0.117 0.906958    
## Male1        0.093835   0.099917   0.939 0.347663    
## Immsupp1    -0.021805   0.104509  -0.209 0.834729    
## copd1       -0.238111   0.145005  -1.642 0.100573    
## chf1         0.041066   0.196049   0.209 0.834081    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.8055) family taken to be 1)
## 
##     Null deviance: 315.11  on 285  degrees of freedom
## Residual deviance: 290.28  on 277  degrees of freedom
## AIC: 1776.9
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.806 
##           Std. Err.:  0.174 
## 
##  2 x log-likelihood:  -1756.851

cbind(IRR = exp(coef(MVdays_nb2)), exp(confint(MVdays_nb2)))[2,]

##      IRR    2.5 %   97.5 % 
## 1.045064 1.020284 1.070437

#Standard errors checked -> okay
MVdays_nbq <- glm.nb(days_MV ~ quartile+APACHE+Age+non_white+Male+Immsupp+copd+chf, data=intub)
summary(MVdays_nbq)

## 
## Call:
## glm.nb(formula = days_MV ~ quartile + APACHE + Age + non_white + 
##     Male + Immsupp + copd + chf, data = intub, init.theta = 1.828470791, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.8183  -0.9216  -0.4349   0.2329   4.2200  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.212944   0.262174   8.441  < 2e-16 ***
## quartile2    0.020217   0.167878   0.120  0.90415    
## quartile3    0.434000   0.161212   2.692  0.00710 ** 
## quartile4    0.449699   0.157266   2.859  0.00424 ** 
## APACHE       0.001820   0.006392   0.285  0.77586    
## Age         -0.008330   0.003165  -2.632  0.00849 ** 
## non_white1   0.015970   0.137671   0.116  0.90765    
## Male1        0.096665   0.100921   0.958  0.33815    
## Immsupp1    -0.027713   0.104094  -0.266  0.79006    
## copd1       -0.275100   0.145398  -1.892  0.05848 .  
## chf1         0.029604   0.195470   0.151  0.87962    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.8285) family taken to be 1)
## 
##     Null deviance: 318.28  on 285  degrees of freedom
## Residual deviance: 289.70  on 275  degrees of freedom
## AIC: 1777.4
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.828 
##           Std. Err.:  0.176 
## 
##  2 x log-likelihood:  -1753.373

#cbind(IRR = exp(coef(MVdays_nbq)), exp(confint(MVdays_nbq)))[2:4,]

#Base = 1
contrasts(intub$quartile) <-contr.treatment(4,base = 1)
MVdays_nb2_q1 <- glm.nb(days_MV ~ quartile+APACHE+Age+non_white+Male+Immsupp+copd+chf, data=intub)
summary(MVdays_nb2_q1)

## 
## Call:
## glm.nb(formula = days_MV ~ quartile + APACHE + Age + non_white + 
##     Male + Immsupp + copd + chf, data = intub, init.theta = 1.828470791, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.8183  -0.9216  -0.4349   0.2329   4.2200  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.212944   0.262174   8.441  < 2e-16 ***
## quartile2    0.020217   0.167878   0.120  0.90415    
## quartile3    0.434000   0.161212   2.692  0.00710 ** 
## quartile4    0.449699   0.157266   2.859  0.00424 ** 
## APACHE       0.001820   0.006392   0.285  0.77586    
## Age         -0.008330   0.003165  -2.632  0.00849 ** 
## non_white1   0.015970   0.137671   0.116  0.90765    
## Male1        0.096665   0.100921   0.958  0.33815    
## Immsupp1    -0.027713   0.104094  -0.266  0.79006    
## copd1       -0.275100   0.145398  -1.892  0.05848 .  
## chf1         0.029604   0.195470   0.151  0.87962    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.8285) family taken to be 1)
## 
##     Null deviance: 318.28  on 285  degrees of freedom
## Residual deviance: 289.70  on 275  degrees of freedom
## AIC: 1777.4
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.828 
##           Std. Err.:  0.176 
## 
##  2 x log-likelihood:  -1753.373

cbind(IRR = exp(coef(MVdays_nb2_q1)),
      exp(confint(MVdays_nb2_q1)),
      p = coef(summary(MVdays_nb2_q1))[,4])

##                   IRR     2.5 %     97.5 %            p
## (Intercept) 9.1425923 5.3009692 15.8962705 3.153570e-17
## quartile2   1.0204223 0.7272564  1.4256791 9.041472e-01
## quartile3   1.5434189 1.1091160  2.1368015 7.100218e-03
## quartile4   1.5678399 1.1401674  2.1435081 4.243256e-03
## APACHE      1.0018216 0.9886235  1.0151979 7.758557e-01
## Age         0.9917047 0.9851211  0.9982833 8.489914e-03
## non_white1  1.0160986 0.7766884  1.3417464 9.076497e-01
## Male1       1.1014916 0.8993489  1.3480047 3.381473e-01
## Immsupp1    0.9726676 0.7968709  1.1890068 7.900622e-01
## copd1       0.7594960 0.5735750  1.0140200 5.848437e-02
## chf1        1.0300465 0.7087438  1.5327841 8.796204e-01

#Base = 2
contrasts(intub$quartile) <-contr.treatment(4,base = 2)
MVdays_nb2_q2 <- glm.nb(days_MV ~ quartile+APACHE+Age+non_white+Male+Immsupp+copd+chf, data=intub)
summary(MVdays_nb2_q2)

## 
## Call:
## glm.nb(formula = days_MV ~ quartile + APACHE + Age + non_white + 
##     Male + Immsupp + copd + chf, data = intub, init.theta = 1.828470791, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.8183  -0.9216  -0.4349   0.2329   4.2200  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.233161   0.274568   8.133 4.18e-16 ***
## quartile1   -0.020217   0.167878  -0.120  0.90415    
## quartile3    0.413783   0.137860   3.001  0.00269 ** 
## quartile4    0.429482   0.134143   3.202  0.00137 ** 
## APACHE       0.001820   0.006392   0.285  0.77586    
## Age         -0.008330   0.003165  -2.632  0.00849 ** 
## non_white1   0.015970   0.137671   0.116  0.90765    
## Male1        0.096665   0.100921   0.958  0.33815    
## Immsupp1    -0.027713   0.104094  -0.266  0.79006    
## copd1       -0.275100   0.145398  -1.892  0.05848 .  
## chf1         0.029604   0.195470   0.151  0.87962    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.8285) family taken to be 1)
## 
##     Null deviance: 318.28  on 285  degrees of freedom
## Residual deviance: 289.70  on 275  degrees of freedom
## AIC: 1777.4
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.828 
##           Std. Err.:  0.176 
## 
##  2 x log-likelihood:  -1753.373

cbind(IRR = exp(coef(MVdays_nb2_q2)),
      exp(confint(MVdays_nb2_q2)),
      p = coef(summary(MVdays_nb2_q2))[,4])

##                   IRR     2.5 %     97.5 %            p
## (Intercept) 9.3293051 5.3385349 16.4324192 4.175382e-16
## quartile1   0.9799864 0.7014201  1.3750308 9.041472e-01
## quartile3   1.5125295 1.1496822  1.9884426 2.686837e-03
## quartile4   1.5364618 1.1785024  2.0001430 1.366298e-03
## APACHE      1.0018216 0.9886235  1.0151979 7.758557e-01
## Age         0.9917047 0.9851211  0.9982833 8.489914e-03
## non_white1  1.0160986 0.7766884  1.3417464 9.076497e-01
## Male1       1.1014916 0.8993489  1.3480047 3.381473e-01
## Immsupp1    0.9726676 0.7968709  1.1890068 7.900622e-01
## copd1       0.7594960 0.5735750  1.0140200 5.848437e-02
## chf1        1.0300465 0.7087438  1.5327841 8.796204e-01

#Base = 3
contrasts(intub$quartile) <-contr.treatment(4,base = 3)
MVdays_nb2_q3 <- glm.nb(days_MV ~ quartile+APACHE+Age+non_white+Male+Immsupp+copd+chf, data=intub)
summary(MVdays_nb2_q3)

## 
## Call:
## glm.nb(formula = days_MV ~ quartile + APACHE + Age + non_white + 
##     Male + Immsupp + copd + chf, data = intub, init.theta = 1.828470791, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.8183  -0.9216  -0.4349   0.2329   4.2200  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.646944   0.269969   9.805  < 2e-16 ***
## quartile1   -0.434000   0.161212  -2.692  0.00710 ** 
## quartile2   -0.413783   0.137860  -3.001  0.00269 ** 
## quartile4    0.015699   0.123956   0.127  0.89922    
## APACHE       0.001820   0.006392   0.285  0.77586    
## Age         -0.008330   0.003165  -2.632  0.00849 ** 
## non_white1   0.015970   0.137671   0.116  0.90765    
## Male1        0.096665   0.100921   0.958  0.33815    
## Immsupp1    -0.027713   0.104094  -0.266  0.79006    
## copd1       -0.275100   0.145398  -1.892  0.05848 .  
## chf1         0.029604   0.195470   0.151  0.87962    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.8285) family taken to be 1)
## 
##     Null deviance: 318.28  on 285  degrees of freedom
## Residual deviance: 289.70  on 275  degrees of freedom
## AIC: 1777.4
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.828 
##           Std. Err.:  0.176 
## 
##  2 x log-likelihood:  -1753.373

cbind(IRR = exp(coef(MVdays_nb2_q3)),
      exp(confint(MVdays_nb2_q3)),
      p = coef(summary(MVdays_nb2_q3))[,4])

##                    IRR     2.5 %     97.5 %            p
## (Intercept) 14.1108494 7.9937450 25.0773939 1.075494e-22
## quartile1    0.6479123 0.4679892  0.9016189 7.100218e-03
## quartile2    0.6611441 0.5029062  0.8698056 2.686837e-03
## quartile4    1.0158227 0.7939350  1.2988435 8.992187e-01
## APACHE       1.0018216 0.9886235  1.0151979 7.758557e-01
## Age          0.9917047 0.9851211  0.9982833 8.489914e-03
## non_white1   1.0160986 0.7766884  1.3417464 9.076497e-01
## Male1        1.1014916 0.8993489  1.3480047 3.381473e-01
## Immsupp1     0.9726676 0.7968709  1.1890068 7.900622e-01
## copd1        0.7594960 0.5735750  1.0140200 5.848437e-02
## chf1         1.0300465 0.7087438  1.5327841 8.796204e-01

#Base = 4
contrasts(intub$quartile) <-contr.treatment(4,base = 4)
MVdays_nb2_q4 <- glm.nb(days_MV ~ quartile+APACHE+Age+non_white+Male+Immsupp+copd+chf, data=intub)
summary(MVdays_nb2_q4)

## 
## Call:
## glm.nb(formula = days_MV ~ quartile + APACHE + Age + non_white + 
##     Male + Immsupp + copd + chf, data = intub, init.theta = 1.828470791, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.8183  -0.9216  -0.4349   0.2329   4.2200  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.662643   0.263364  10.110  < 2e-16 ***
## quartile1   -0.449699   0.157266  -2.859  0.00424 ** 
## quartile2   -0.429482   0.134143  -3.202  0.00137 ** 
## quartile3   -0.015699   0.123956  -0.127  0.89922    
## APACHE       0.001820   0.006392   0.285  0.77586    
## Age         -0.008330   0.003165  -2.632  0.00849 ** 
## non_white1   0.015970   0.137671   0.116  0.90765    
## Male1        0.096665   0.100921   0.958  0.33815    
## Immsupp1    -0.027713   0.104094  -0.266  0.79006    
## copd1       -0.275100   0.145398  -1.892  0.05848 .  
## chf1         0.029604   0.195470   0.151  0.87962    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.8285) family taken to be 1)
## 
##     Null deviance: 318.28  on 285  degrees of freedom
## Residual deviance: 289.70  on 275  degrees of freedom
## AIC: 1777.4
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.828 
##           Std. Err.:  0.176 
## 
##  2 x log-likelihood:  -1753.373

cbind(IRR = exp(coef(MVdays_nb2_q4)),
      exp(confint(MVdays_nb2_q4)),
      p = coef(summary(MVdays_nb2_q4))[,4])

##                    IRR     2.5 %     97.5 %            p
## (Intercept) 14.3341208 8.3166527 24.8848017 4.981497e-24
## quartile1    0.6378202 0.4665249  0.8770642 4.243256e-03
## quartile2    0.6508460 0.4999643  0.8485345 1.366298e-03
## quartile3    0.9844238 0.7699157  1.2595489 8.992187e-01
## APACHE       1.0018216 0.9886235  1.0151979 7.758557e-01
## Age          0.9917047 0.9851211  0.9982833 8.489914e-03
## non_white1   1.0160986 0.7766884  1.3417464 9.076497e-01
## Male1        1.1014916 0.8993489  1.3480047 3.381473e-01
## Immsupp1     0.9726676 0.7968709  1.1890068 7.900622e-01
## copd1        0.7594960 0.5735750  1.0140200 5.848437e-02
## chf1         1.0300465 0.7087438  1.5327841 8.796204e-01

#Reset base level
contrasts(intub$quartile) <-contr.treatment(4,base = 1)

#ICU days (mean x IRR)
incr_days = (1.058129 * mean(imgs$ICU_days)) - mean(imgs$ICU_days)
incr_days

## [1] 0.4002597

median(imgs$ICU_days)

## [1] 5

#Days MV
more_days = (1.045064 * mean(intub$days_MV)) - mean(intub$days_MV)
more_days

## [1] 0.3691782

median(intub$days_MV)

## [1] 5

library(reshape2)
library(RColorBrewer)
#Eliminate non-intubated from the MV data and transform for plotting
meltICU <- imgs %>% dplyr::select(Subject_ID, ICU_days, quartile) %>% melt(id.vars="quartile", measure.vars="ICU_days")
meltMV <- imgs %>% filter(days_MV != 0) %>% dplyr::select(Subject_ID, days_MV, quartile) %>% melt(id.vars="quartile", measure.vars="days_MV")
meltplot <- rbind(meltICU, meltMV)

ICUdata <- imgs %>% dplyr::select(Subject_ID, ICU_days, quartile) %>% rename(value = ICU_days) %>% mutate(variable = "ICU_days")
MVdata <- intub %>% dplyr::select(Subject_ID, days_MV, quartile) %>% rename(value = days_MV) %>% mutate(variable = "days_MV")
meltplot <- rbind(ICUdata, MVdata)

#Boxplot square root scale
ggplot(meltplot) + geom_boxplot(aes(quartile, sqrt(value), color=variable)) +   
  scale_x_discrete(name="Chest Radiograph Score Quartile") + 
  scale_y_continuous(name="(Square root of) Days") +
  theme_minimal() + theme(plot.title = element_text(face = "bold")) +
  theme(panel.grid.minor = element_blank()) + 
  theme(legend.title = element_blank()) +
  theme(legend.position = "bottom") +
  scale_color_discrete(labels=c("Length of stay", "Ventilator days"))

#Boxplot log scale
ggplot(meltplot) + geom_boxplot(aes(quartile, log(value), color=variable)) +   
  scale_x_discrete(name="Chest Radiograph Score Quartile") + 
  scale_y_continuous(name="(Log of) Days", limits=c(0,5)) +
  theme_minimal() + theme(plot.title = element_text(face = "bold")) +
  theme(panel.grid.minor = element_blank()) + 
  theme(legend.title = element_blank()) +
  theme(legend.position = "bottom") +
  scale_color_discrete(labels=c("Length of stay", "Ventilator days"))

box_comparisons <- list(c("1","2"), c("2","3"), c("3","4"), c("1","3"), c("1","4"),  c("2","4"))

symnum.args <- c("***"=0.05,ns=1)

var_names <- c('ICU_days' = "Length of ICU stay", 'days_MV' = "Ventilator days")

my_boxplot <- ggplot(meltplot, aes(quartile, log(value))) + geom_boxplot() + 
  facet_grid(.~ variable, labeller=as_labeller(var_names)) +
  my_stat_compare_means(label = "p.signif", method = "t.test", comparisons = box_comparisons, symnum.args = symnum.args) +
  scale_x_discrete(name="Chest Radiograph Score Quartile") +
  scale_y_continuous(name="(Log) Days", limits=c(-0.25,7.5)) +
  theme_bw() +
  annotate("text", x=2.5, y=-0.2, label="p for trend < .001")

my_boxplot 

wilcox_data <- meltplot %>% filter(variable == "days_MV") %>% filter(quartile == 2 | quartile == 3)
t.test(log(wilcox_data$value) ~ wilcox_data$quartile)

## 
##  Welch Two Sample t-test
## 
## data:  log(wilcox_data$value) by wilcox_data$quartile
## t = -2.0012, df = 148.09, p-value = 0.0472
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.56153943 -0.00353892
## sample estimates:
## mean in group 2 mean in group 3 
##        1.512006        1.794545

# jonkheere Terpstra test
ICUdata$quartile <- ordered(ICUdata$quartile, levels=c("1", "2", "3", "4"))
jonckheere.test(ICUdata$value, ICUdata$quartile, alternative="two.sided")

## Warning in jonckheere.test(ICUdata$value, ICUdata$quartile, alternative = "two.sided"): Sample size > 100 or data with ties 
##  p-value based on normal approximation. Specify nperm for permutation p-value

## 
##  Jonckheere-Terpstra test
## 
## data:  
## JT = 74022, p-value = 1.16e-12
## alternative hypothesis: two.sided

MVdata$quartile <- ordered(MVdata$quartile, levels=c("1", "2", "3", "4"))
jonckheere.test(MVdata$value, MVdata$quartile, alternative="two.sided")

## Warning in jonckheere.test(MVdata$value, MVdata$quartile, alternative = "two.sided"): Sample size > 100 or data with ties 
##  p-value based on normal approximation. Specify nperm for permutation p-value

## 
##  Jonckheere-Terpstra test
## 
## data:  
## JT = 18238, p-value = 3.097e-05
## alternative hypothesis: two.sided

Miscellaneous

shapiro.test(imgs$survival)

## 
##  Shapiro-Wilk normality test
## 
## data:  imgs$survival
## W = 0.79791, p-value < 2.2e-16

imgs %>% summarize(mean = mean(survival), median = median(survival), min = min(survival), max = max(survival), count = length(survival))

## # A tibble: 1 x 5
##    mean median   min   max count
##   <dbl>  <dbl> <dbl> <dbl> <int>
## 1  932.   364.     0  3664   560

imgs %>% mutate(surv_quart = ntile(survival, 4)) %>% group_by(surv_quart) %>% summarize(mean = mean(survival), median = median(survival), min = min(survival), max = max(survival), count = length(survival))

## # A tibble: 4 x 6
##   surv_quart    mean median   min   max count
##        <int>   <dbl>  <dbl> <dbl> <dbl> <int>
## 1          1    8.89     8      0    22   140
## 2          2  127.     104     22   362   140
## 3          3  951.     912.   365  1788   140
## 4          4 2643.    2388   1796  3664   140

imgs %>% group_by(Final_Dx) %>% summarize(count = length(Total))

## # A tibble: 4 x 2
##   Final_Dx    count
##   <chr>       <int>
## 1 ARDS           26
## 2 Sepsis        292
## 3 Sepsis/ARDS    97
## 4 SIRS          145

imgs %>% ggplot(aes(Total, fill = quartile)) + geom_histogram(bins = 15)

imgs %>% summarize( median = median(Total), mean = mean(Total))

## # A tibble: 1 x 2
##   median  mean
##    <dbl> <dbl>
## 1      7  6.61

imgs %>% group_by(quartile) %>% summarize(min = min(Total), max = max(Total))

## # A tibble: 4 x 3
##   quartile   min   max
##   <fct>    <dbl> <dbl>
## 1 1            0     2
## 2 2            2     7
## 3 3            7    11
## 4 4           11    16

imgs %>% summarize(med_age = median(Age), med_apa = median(APACHE))

## # A tibble: 1 x 2
##   med_age med_apa
##     <dbl>   <dbl>
## 1      60      25

imgs %>% group_by(Male) %>% summarize(count = length(Total))

## # A tibble: 2 x 2
##   Male  count
##   <fct> <int>
## 1 0       260
## 2 1       300

imgs %>% group_by(non_white) %>% summarize(count = length(Total))

## # A tibble: 2 x 2
##   non_white count
##   <fct>     <int>
## 1 0           437
## 2 1           123

imgs %>% group_by(Immsupp) %>% summarize(count = length(Total))

## # A tibble: 2 x 2
##   Immsupp count
##   <fct>   <int>
## 1 0         355
## 2 1         205

imgs %>% mutate(age_qt = ntile(Age,4)) %>% group_by(age_qt) %>% summarize(lower = min(Age), upper = max(Age))

## # A tibble: 4 x 3
##   age_qt lower upper
##    <int> <dbl> <dbl>
## 1      1    19    47
## 2      2    47    60
## 3      3    60    70
## 4      4    70    94

imgs %>% mutate(apa_qt = ntile(APACHE,4)) %>% group_by(apa_qt) %>% summarize(lower = min(APACHE), upper = max(APACHE))

## # A tibble: 4 x 3
##   apa_qt lower upper
##    <int> <dbl> <dbl>
## 1      1     3    20
## 2      2    20    25
## 3      3    25    31
## 4      4    31    54

imgs %>% filter(!is.na(date_death)) %>% summarize(count = length(Total)) #total death events in cohort (Cox event rate)

## # A tibble: 1 x 1
##   count
##   <int>
## 1   369

imgs %>% filter(DC_coded == 0 | DC_coded == 1) %>% summarize(count = length(Total)) #in hospital deaths - entire cohort

## # A tibble: 1 x 1
##   count
##   <int>
## 1   159

imgs_non %>% filter(DC_coded == 0 | DC_coded == 1) %>% summarize(count = length(Total)) # in hospital deaths for ARDS only

## # A tibble: 1 x 1
##   count
##   <int>
## 1    93

key <- read_excel("~/Documents/ROCIanalysis/ROCI_101918.xlsx", sheet = "Key", range = cell_cols(1:2))
subjects <- key %>% filter(Subject_ID %in% imgs$Subject_ID)
subjects <- imgs %>% dplyr::select(Subject_ID, cxr_date, Hosp_dc) %>% right_join(subjects)
write.csv(subjects, file="subjects.csv", row.names = FALSE)

#CXR scores are non-normally distributed and ordinal, will use a spearman correlation.

imgs %>% dplyr::select(Total, APACHE) %>% cor(., use = "pairwise.complete.obs", method="spearman")

##            Total    APACHE
## Total  1.0000000 0.1764567
## APACHE 0.1764567 1.0000000

cor.test(imgs$Total, imgs$APACHE, method="spearman", alternative = "t")

## Warning in cor.test.default(imgs$Total, imgs$APACHE, method = "spearman", :
## Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  imgs$Total and imgs$APACHE
## S = 24104000, p-value = 2.676e-05
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.1764567

imgs %>% ggplot(aes(Subject_ID, Total, color=hosp_death)) + geom_jitter() + scale_color_manual(values=c("black", "red"))

imgs %>% ggplot(aes(quartile, Total, color=hosp_death)) + geom_jitter() + scale_color_manual(values=c("black", "red"))

imgs %>% ggplot(aes(Total, fill = hosp_death)) + geom_histogram() + scale_fill_manual(values=c("black", "red"))